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Abstract 

Shock  tunnels  and  expansion  tubes  are  our  tools  of  choice  for  producing  high-enthalpy  flows  in 
aerothermodynamic  studies  relevant  to  hypersonic  flight.  The  detailed  flow  in  these  machines  turns  out 
to  be  quite  complex  and  computer  models  are  needed  to  completely  and  accurately  describe  the  experi¬ 
mental  flow  properties  and  to  diagnose  the  gas  dynamic  behaviour  of  the  machine  when  developing  new 
operating  conditions.  We  describe  a  couple  of  the  flow  simulation  codes  that  we  have  written  to  assist 
us  in  our  experimental  work.  The  first  is  a  quasi-one-dimensional  flow  code  that  is  capable  of  modelling 
entire  free-piston  driven  facilities,  albeit  with  some  significant  limitations.  The  second  is  a  finite-volume 
flow  code  that  can  more  accurately  capture  the  strong  viscous  and  thermochemical  interactions  that  affect 
strongly  expanding  flows  that  are  an  operating  characteristic  of  expansion  tube  facilities.  Some  applica¬ 
tion  studies  are  then  described. 


1  INTRODUCTION 

High  flow  speeds,  in  the  range  of  4  to  10  km/s,  are  required  for  the  ground-based  study  of  supersonic  combus¬ 
tion  and  the  aerothermodynamics  of  planetary  entry.  The  operation  of  a  continuous  wind  tunnel  producing 
the  requisite  flow  speeds  is  problematic.  First,  it  is  difficult  to  design  a  facility  that  would  bear  the  sustained 
heating  rates  to  the  tube  walls  for  long  periods  and,  second,  the  energy  requirements  for  continuous  opera¬ 
tion  are  prohibitive.  Impulse  facilities,  such  as  shock  tunnels  and  expansion  tubes,  arc  thus  the  machines  of 
choice  for  testing  at  hypersonic  speeds. 

The  hypersonics  research  group  at  The  University  of  Queensland  has  a  number  of  free-piston  driven 
shock  tunnels  and  expansion  tubes  designed  to  produce  high-enthalpy,  short-duration  test  flows.  These 
machines  have  been  constructed  and  operated  over  the  past  three  decades  and  there  is  now  a  degree  of 
confidence  in  designing  new  machines  and  predicting  their  behaviour.  Although  the  gas  dynamic  processes 
in  the  machines  are  dominated  by  unsteady  one-dimensional  wave  interactions  which  can  be  sketched  as  a 
wave  diagram,  there  are  significant  secondary  processes  driven  by  viscous  effects  and  finite -rate  chemical 
kinetics.  To  cope  with  increasing  model  complexity,  a  number  of  flow  simulation  tools  have  been  developed 
over  the  past  twenty  years  to 

•  assist  the  design  of  new  machine  components  such  as  hypersonic  nozzles, 

•  provide  the  answers  to  “What  if...”  questions  in  the  design  of  new  operating  conditions,  and  to 

•  be  a  key  component  in  fully  specifying  each  facility’s  test  flow  conditions. 
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ORGANIZATION 


This  paper  will  be  limited  to  describing  our  efforts  in  modelling  the  Drummond  shock  tunnel  and  the  X2 
expansion  tube,  as  well  as  the  computer  codes  that  we  use  for  that  modelling.  These  particular  facilities  have 
been  selected  to  demonstrate  different  capabilities  of  the  simulation  tools.  The  Drummond  shock  tunnel 
is  a  good  starting  example  because  relatively  simple  computational  modelling  can  do  a  reasonable  job  of 
estimating  its  overall  performance.  Also,  there  arc  complex  fluid  dynamic  interactions  associated  with  the 
shock  reflection  process  that  can  be  explored  in  the  absence  of  high-temperature  thermochemical  effects.  The 
X2  facility  is  a  good  example  of  a  high  performance  machine,  and  comes  with  the  associated  complications 
of  a  free-piston  driver  and  significant  thermochemical  nonequilibrium. 


2  THE  FLOW  FACILITIES 

2.1  Small  Shock  Tunnel 

The  Small  Shock  Tunnel  (SST)  facility  (Figure  1)  operated  by  the  hypersonics  group  at  UQ  is  an  old  shock 
tube  from  the  Australian  Defence  Science  and  Technology  Organisation  (DSTO)  [1]  that  has  been  redevel¬ 
oped  into  a  reflected-shock  tunnel.  It  is  also  known  as  the  Drummond  Tunnel  after  the  DSTO  researcher 
responsible  for  its  original  construction.  The  tunnel  has  been  used  as  (i)  a  facility  for  testing  new  ideas  for 
hypersonic  nozzles  [2]  and  diagnosing  tunnel  operation  problems  [3] ;  (ii)  a  relatively  low  enthalpy  gas  dy¬ 
namics  facility  for  both  teaching  and  research;  and  (iii)  a  test  facility  for  the  development  of  laser-based  flow 
diagnostics  [4]. 


Figure  1:  Layout  of  the  small  reflected-shock  tube  known  as  the  Drummond  tunnel  [2]. 

The  motivation  for  using  the  SST  facility  as  a  prototyping  facility  arose  from  an  anomoly  seen  in  the 
flows  produced  by  large  free-piston  driven  shock  tubes.  Calibrations  of  Mach  8  axisymmetric  nozzles  on 
both  the  T4  facility  at  UQ  [5]  and  the  HEG  shock  tunnel  in  Gottingen  [6]  showed  significant  disturbances 
in  Pitot  pressure  near  the  centreline  of  the  test  flow.  These  disturbances  occur  at  the  start  of  the  flow  period 
and,  for  high  enthalpies,  make  much  of  the  test  flow  unusable.  The  SST  facility  was  used,  together  with  the 
mbcns2  code,  to  identify  the  flow  disturbance  mechanism. 

The  SST  initally  operated  in  a  non-reflecting  configuration  in  which  the  test  gas  was  accelerated  by  an 
incident  shock  and  then  flowed  over  the  model  under  test.  This  produced  a  high  speed  but  relatively  low 
Mach  number  ( M  <  2)  flow  in  a  tube  with  limited  optical  access.  In  order  to  produce  flows  with  higher 
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Mach  numbers,  the  facility  was  modified  to  include  a  converging-diverging  nozzle  and  a  test  section  with 
reasonably  large  windows.  With  these  additions,  the  facility  is  now  operated  in  a  reflected-shock  mode  where 
the  nozzle  effectively  closes  off  the  downstream  end  of  the  shock  tube.  As  shown  in  Figure  2,  the  test  gas 
that  is  processed  by  the  incident  shock  is  brought  to  rest  upstream  of  the  nozzle  throat  via  a  reflected  shock. 
This  hot,  nearly  stagnant  gas  is  then  expanded  through  the  converging-diverging  nozzle  to  a  Mach  number 
of  approximately  4  or  7,  depending  on  the  nozzle  contour.  The  static  temperature  in  the  test  flow  is  quite  low 
for  either  nozzle  so,  for  generating  flows  at  the  high-end  of  the  enthalpy  range,  we  typically  use  a  facility 
with  a  free-piston  driver. 


Tailored 

Time  / 


Figure  2:  Wave  diagram  and  internal  view  of  the  Drummond  tunnel  [2]. 


2.2  X2  Expansion  Tube 

There  are  three  expansion  tube  facilities  at  the  University  of  Queensland.  In  order  of  their  commissioning, 
and  physical  size,  they  are:  XI,  X2  and  X3.  The  following  discussion  will  focus  on  the  specifics  of  the  X2 
facility,  though  the  principles  of  operation  are  the  same  in  all  three  facilities.  Much  of  their  recent  work  has 
been  directed  toward  the  study  of  the  aerothermodynamics  of  naturally  radiating  flows. 

In  expansion  tube  mode,  experiments  are  performed  on  subscale  models  of  aeroshell  configurations. 
Because  the  models  in  the  experiment  are  scaled,  the  flow  conditions  also  require  appropriate  scaling  to 
maintain  similarity  between  the  experiment  and  the  true  flight  conditions.  Typically,  the  product  of  density 
and  length  is  kept  the  same  as  for  the  full  scale  vehicle  and,  although  this  provides  similarity  for  some 
aerothermodynamic  processes,  the  radiative  energy  exchange  scales  differently.  When  it  is  necessary  to 
study  a  radiating  flow  field  at  exact  flight  conditions,  the  expansion  tube  may  be  operated  in  nonrcflcctcd 
shock  tube  mode.  When  using  nonreflected  shock  tube  mode,  there  is  no  subscale  model  placed  in  the  test 
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flow.  Instead,  the  object  of  observation  is  the  relaxing  flow  immediately  behind  the  shock  as  it  propagates 
into  the  test  gas.  It  is  helpful  to  think  of  this  type  of  experiment  as  observing  the  bow  shock  that  is  driven  by 
the  aeroshell.  There  is  no  actual  aeroshell  driving  the  shock,  rather  just  high  pressure  gas  generated  in  the 
shock  tube.  The  two  types  of  experiment,  subscale  testing  in  expansion  tube  mode  and  direct  observation  in 
nonreflected  shock  tube  mode,  are  shown  in  Figure  3. 


radiating  sublayer 
at  true-flight  conditions 


Figure  3:  Schematic  representation  of  types  of  experiments  performed  in  X2  for  the  study  of  radiating  blunt 
body  flows. 


Figure  4  shows  a  schematic  of  the  X2  facility  and  a  wave  diagram  of  the  flow  process.  The  representation 
in  Figure  4  corresponds  to  use  of  the  facility  in  expansion  tunnel  mode.  The  test  gas  which  eventually  passes 
over  the  model  (or  in  this  diagram,  the  Pitot  rake),  initially  resides  in  the  shock  tube.  The  gas  initially  at  rest 
in  the  acceleration  tube  is  called  the  accelerator  gas.  In  the  diagram,  t  =  0  corresponds  to  the  rupture  of  the 
primary  diaphragm;  at  this  point  the  35  kg  piston  is  nearing  the  end  of  its  stroke.  The  primary  diaphragm 
is  made  of  mild  steel.  Its  thickness  is  varied  according  to  the  desired  rupture  pressure.  The  rupture  of 
the  diaphragm  produces  a  shock,  labelled  the  primary  shock.  This  shock  travels  through  the  shock  tube 
compressing  and  accelerating  the  test  gas.  The  secondary  diaphragm  is  only  very  light  and  so  ruptures  when 
the  primary  shock  arrives.  The  secondary  shock  propagates  into  the  acceleration  tube,  its  faster  speed  a 
consequence  of  the  lower  density  medium.  The  compressed  test  gas  now  expands,  in  an  unsteady  manner, 
into  the  acceleration  tube.  This  expansion  process  drives  the  low  pressure  accelerator  gas  before  it.  This 
interface  between  test  gas  and  accelerator  gas  is  labelled  as  the  contact  surface.  The  test  gas  expansion 
is  controlled  by  the  pressure  of  the  accelerator  gas;  a  lower  pressure  accelerator  gas  results  in  a  stronger 
expansion.  The  steady  test  time,  as  indicated  in  Figure  4,  is  limited  by  the  arrival  of  a  u  +  a  wave  from  the 
centred  expansion.  For  further  details  about  the  principles  of  expansion  tubes  see  References  [7,  8]. 

Numerical  simulations  aid  the  expansion  tube  testing  by  estimating  a  full  set  of  flow  properties  that  are 
not  directly  measurable.  For  example,  numerical  simulations  help  to  establish  the  set  of  flow  conditions  in 
the  test  gas  when  operated  in  expansion  tunnel  mode.  This  is  important  because  these  are  the  flow  conditions 
which  form  the  free  stream  flow  for  the  subscale  model;  subsequent  analysis  of  the  experiment  depends 
on  reliable  estimates  of  the  free  stream  flow  conditions.  An  example  simulation  for  a  condition  in  the  X2 
facility,  operated  in  expansion  tube  mode  is  presented  in  Section  4.4.  This  simulation  targets  the  25  MJ/kg 
CO2-N2  operating  condition  developed  in  Reference  [9]  to  obtain  spectrally  resolved  radiation  intensity 
measurements  of  the  recompression  shock  formed  over  a  finite  cylinder.  The  model  and  optics  layout  in  the 
test-section  is  illustrated  in  Figure  5. 
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Driver  gas  unsteady  expansion  Test  gas  unsteady  expansion 


Figure  4:  A  schematic  of  the  X2  facility  configured  for  expansion  tunnel  operation.  A  wave  diagram  (x  —  t 
diagram)  shows  the  flow  process. 


Figure  5:  Schematic  of  bluntbody  shock  layer  spectroscopy  experiment  with  25mm  OD  by  102nmi  cylinder. 
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The  X2  facility  is  operated  as  a  conventional  shock  tube  when  used  in  nonreflected  shock  tube  mode. 
Figure  6  shows  a  schematic  of  the  facility  and  a  wave  diagram  for  the  flow  process.  Prior  to  t  =  0,  the 
piston  compresses  the  driver  gas  which  is  often  a  mixture  of  helium  and  argon.  The  rupture  of  the  diaphragm 
corresponds  to  t  =  0  in  Figure  6.  After  diaphragm  rupture,  the  shock  then  propagates  into  the  shock  tube 
which  contains  the  test  gas.  The  experiment  involves  making  optical  measurements  of  the  radiating  flow 
as  the  shock  emerges  from  the  shock  tube  and  into  the  test  section.  The  section  of  the  flow  labelled  test 
time  in  Figure  6  distinguishes  this  mode  of  operation  from  the  expansion  tunnel  mode:  for  the  nonreflected 
shock  tube  mode,  the  gas  immediately  behind  the  primary  shock  is  of  interest  whereas  for  the  expansion 
tunnel  mode,  the  test  gas  has  been  previously  “gathered  up”  in  the  shock  tube  and  allowed  to  expand  in  the 
accelerator  tube.  In  the  experiments  using  nonreflected  shock  tube  mode  only  a  small  portion  of  the  test  gas 
is  of  interest,  namely,  the  nonequilibrium  radiating  gas  immediately  behind  the  shock. 
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Figure  6:  A  schematic  of  the  X2  facility  configured  for  nonreflected  shock  tube  operation.  A  wave  diagram 
(x  —  t  diagram)  shows  the  flow  process.  This  figure  is  provided  courtesy  of  Carolyn  Jacobs. 
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3  THE  SIMULATION  CODES 

Because  the  experimental  facilities  are  complex  multiscale  machines,  we  have  a  range  of  codes,  each  with 
quite  different  capability  and  computational  requirements.  The  codes  that  will  be  discussed  here  are: 

•  Lld2:  a  simulation  code  for  time-dependent  quasi-one-dimensional  flows,  used  to  simulate  entire 
facilities  with  moderate  level  of  detail,  and 

•  mbcns2:  a  time-accurate  simulation  code  for  two-dimensional,  axisymmetric  or  planar  flows,  used 
when  multidimensional  effects  cannot  be  separated  from  the  transient  flow. 

3.1  One-dimensional  Flow  Code,  Lld2 

In  order  to  estimate  the  performance  of  a  free -piston  driven  impulse  facility,  one  must  consider  the  both  the 
dynamics  of  the  piston  and  gases,  and  the  viscous  effects  (including  heat  transfer)  simultaneously.  Models 
which  omit  these  effects  require  a  number  of  facility-specific  fudge  factors  which  can  be  obtained  accurately 
only  after  the  construction  and  operation  of  the  facility.  This  section  describes  the  numerical  modelling 
behind  the  computer  code  Lld2,  which  is  capable  of  simulating  the  (gas-dynamic)  operation  of  a  free- 
piston  driven  facility  during  the  design  process.  It  is  closely  related  to  other  light-gas  gun  codes  (see  e.g. 
[10],  [1 1],  [12],  [13],  [14])  and  borrows  a  number  of  ideas  from  some  of  them.  The  principal  features  of  the 
code  are: 

•  Quasi-one-dimensional  formulation  for  the  gas-dynamics.  There  is  only  one  spatial  coordinate  but 
gradual  variation  of  duct  area  is  allowed. 

•  The  ability  to  simulate  several  independent  (or  interacting)  slugs  of  gas.  Also,  several  pistons/pro¬ 
jectiles  and  multiple  diaphragms  may  be  included.  Coupling  to  the  gas  dynamics  is  via  the  boundary 
conditions  of  the  gas  slugs. 

•  A  Lagrangian  discretization  of  the  gas  slugs.  This  is  done  by  dividing  each  gas  slug  into  a  set  of 
control-masses  (or  gas  particles)  and  following  the  positions  of  these  particles. 

•  Nominal  second-order  accuracy  in  both  space  and  time  combined  with  a  robust  shock-capturing  scheme. 
The  use  of  a  shock  capturing  scheme  means  that  the  same  set  of  equations  is  used  to  compute  the  mo¬ 
tion  of  the  gas  whether  a  shock  is  present  or  not.  This  simplifies  the  code  (as  shocks  do  not  need  to  be 
explicitly  identified  or  tracked)  and  is  especially  important  in  situations  where  shocks  may  form  from 
the  merger  of  finite  compression  waves  and  where  multiple  shocks  and  contact  surfaces  interact  in  a 
complicated  manner.  It  also  results  in  a  smearing  of  the  shocks  over  a  couple  of  computational  cells. 
However,  in  practice,  this  is  not  a  problem  as  any  smeared  shocks  can  be  sharpened  by  increasing  the 
grid  resolution. 

•  Viscous  effects  are  included  using  the  standard  engineering  correlations  for  friction  and  heat  transfer 
in  pipe  flow.  Although  these  correlations  are  generally  derived  for  steady  incompressible  flow,  they 
seem  to  perform  adequately  in  the  simulations  where  the  flows  are  predominantly  unsteady  and  are 
very  compressible. 

The  general  procedure  for  modelling  a  specific  facility  (or  system)  is  to  divide  the  facility  into  its  compo¬ 
nent  parts  such  as  the  tube,  pistons,  diaphragms  and  volumes  of  gas  (i.e.  gas  slugs).  The  description  of  each 
component  is  formulated  separately  and  components  allowed  to  interact  through  boundary  conditions.  The 
core  of  LI d2  is  a  time-stepping  loop  which  first  applies  the  specified  boundary  conditions  and  then  advances 
the  state  of  the  entire  system  forward  in  time  by  a  small  increment  (or  time  step).  The  generic  components 
described  in  the  following  sections  include  a  slug  of  compressible  gas,  a  piston  and  a  diaphragm. 
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Figure  7:  A  typical  control-mass  or  Lagrangian  cell. 


3.1.1  Gas  Dynamics 


Each  slug  of  gas  is  treated  in  a  Lagrangian  framework  in  which  the  slug  is  divided  into  a  number  of  control- 
mass  elements  (or  cells)  moving  in  a  variable-area  duct.  Flow  in  one  dimension  only  is  considered  and  any 
area  changes  in  the  tube  area  are  assumed  to  be  gradual.  Although  the  boundary  layer  along  the  tube  wall  is 
not  completely  modelled  in  the  formulation  of  the  gas-dynamic  equations,  some  of  its  effects  are  modelled 
in  the  momentum  equation  by  the  addition  of  a  wall  shear  stress.  These  approximations  are  arguably  the 
most  troublesome  part  of  the  modelling  process  as  they  cannot  be  fixed  later  by  simply  increasing  the  grid 
resolution. 


Figure  7  shows  a  typical  control-mass  cell  (labelled  j)  with  interfaces  (labelled  j  —  \  and  j  +  |)  to 
adjacent  cells.  At  each  interface,  the  Lagrangian  description  equates  the  local  fluid  velocity  to  the  interface 
velocity  as 


dxi±\ 

dt 


~  uiH 


(1) 


where  x  is  the  position  of  the  interface  and  u  is  the  local  gas  velocity  computed  with  a  Riemann  solver  (to 
be  described  later). 


The  average  density  within  the  cell  is  given  by 

_  m7- 

Pi  =  —( - y  * 


(2) 


where  (•  •  •)  represents  a  cell  average,  A  is  the  area  of  the  duct  and  rrij  is  the  (constant)  mass  of  gas  in  cell  j. 

The  rate  of  change  of  momentum  in  the  cell  is  due  to  the  pressure  forces  acting  on  the  cell  interfaces  and 
viscous  forces  acting  at  the  duct  wall.  It  is  given  by 


—  m.nU 

dt 


3  u3 


■P'wall 


(3) 


where  Fwau  is  the  shear  friction  force  at  the  wall  and  Fioss  is  an  effective  body-force  due  to  pipe-fitting 
losses,  for  example.  Evaluation  of  these  loss  terms  will  be  discussed  in  the  viscous  effects  section  3.1.2. 
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The  rate  of  change  of  energy  within  the  cell  is  due  to  the  work  done  at  the  cell  interfaces  plus  the  heat 
transferred  through  the  duct  wall.  It  is  given  by 


dt 


mjEj 


P  i_A 
J  2 


-  p,+iAmui- n+qi 


(4) 


where  E  =  e  +  \u2  is  the  total  specific  energy  of  the  gas  and  q  is  the  rate  of  heat  transfer  into  the  cell. 
Evaluation  of  q  will  appear  later  in  section  3.1.2.  Note  that  there  is  no  shear  stress  term  in  the  total  energy 
equation.  This  is  because  the  energy  removed  from  the  kinetic  energy  component  by  the  wall  shear  stress  is 
deposited  into  the  internal  energy  component  of  the  gas  near  the  wall  and,  since  there  is  no  transfer  of  mass 
from  one  cell  to  the  next,  the  total  energy  of  the  cell  is  unchanged  (by  this  mechanism). 

The  governing  differential  equations  for  the  gas  dynamics  (i.e.  equations  (1),  (3)  and  (4) )  are  completed 
by  specifying  the  thermodynamic  properties  of  the  gas  as  described  in  Section  3.3. 


3.1.2  Viscous  Effects 

The  viscous  shear  force  on  a  gas  cell  is  given  by 

Fwall  =  T0  7T  D  (. Xj+l  -  Xj_l)  ,  (5) 

where  tq  is  the  local  shear  stress  at  the  wall  and  D  is  the  (average)  effective  diameter  of  the  tube.  Assuming 
a  circular  cross-section 

D  =  2{A/tt)1/2  .  (6) 

The  wall  shear  stress  is  obtained  from  the  Darcy  formula  for  for  steady  incompressible  flow  (see  e.g.  [15,  16]) 


TO 


-pfu\u 

8 


(V) 


where  /  is  the  Darcy -Weisbach  friction  factor.  With  minor  changes,  we  follow  [13]  and  use 


/ 

/ 

/ 


64 

A  Re 


,  Re  <  2000  , 


0.032 

A 


Re 

2000 


n  0.3187 


,  2000  <  Re  <  4000 


1 

A 


1.14-2 


l°gto  ( 


21.25  Re~°'9  +  — 


D  J  \ 


1  -2 


,  Re  >  4000 


(8) 


where  the  Re  is  the  local  Reynolds  number  based  on  tube  diameter 


Re  = 


p*  D  |rt| 

p* 


(9) 


and  e  is  the  absolute  wall  roughness.  The  explicit  expression  for  /  for  the  turbulent  regime  is  taken  from 
[17]  and  is  within  1%  of  the  well  known  Colebrook- White  equation.  For  Reynolds  numbers  up  to  105  in 
shock-tube  type  flows,  it  is  reasonable  to  assume  a  smooth  wall  and  use 


/  =  ]-[1.81og10(i?e)-1.5147]-2 


Re  >  4000 


(10) 


The  properties  p*  and  p*  =  pT/T*  are  evaluated  at  the  Eckert  reference  temperature  (see  e.g.  [18]  section 
5.12) 

T*  =  T  +  0.5(TW  —  T)  +  0.22(Taw  —  T)  ,  (11) 
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where  T  is  the  cell-average  temperature,  Tw  is  the  specified  wall  temperature  and 


T„,„  =  AT 


(12) 


is  the  adiabatic  wall  temperature.  Since,  we  are  interested  in  flows  which  may  have  very  high  Mach  numbers, 
a  compressibility  correction  is  applied  via  the  the  compressibility  factor  [19] 


a  =  i  +  ^ ^  n  m2  , 
2 


(13) 


where  M  is  the  local  Mach  number  and  (2  is  the  recovery  factor.  Although  the  compressibility  factor  A  in 
equations  (8)  and  (10)  was  suggested  for  rough  surfaces,  we  have  used  it  to  adjust  the  friction  factor  for 
all  values  of  Re.  For  laminar  flow  (i.e.  Re  <  2000),  the  recovery  factor  is  set  to  (2  =  (Pr)1/2  while,  for 
turbulent  flow,  (2  =  (Pr)1/3. 

Pressure  losses  due  to  sudden  changes  in  tube  cross-section  are  computed  for  each  cell  as 


Floss 


^ Floss  ~r  / 

-l^t  A  {x’+i 


-XJ-V 


(14) 


where  ^ 

APzoss  =  —Kl  -pu\u 


(15) 


and  Lioss  is  the  length  of  tube  over  which  the  pressure  loss  is  distributed.  Values  of  K l  / Lioss  are  stored  along 
with  the  cross-sectional  area  for  the  tube.  For  a  contraction  and  a  diaphragm  station,  we  use  Kl  =  0.25. 


Fleat  transfer  into  a  gas  cell  is  given  by  ([18],  section  5.12) 


q  =  fin  D  (xj+i  -  Xj_i)  (Tw  -  Taw)  ,  (16) 

where  the  heat  transfer  coefficient  is 

h  =  p  Cp  |«|  St  (17) 

and  the  Stanton  number  is  given  by  the  modified  Reynolds  analogy  for  turbulent  flow  in  pipes  ([18],  section 

6.2) 

St  =  {  Pr-2/3  .  (18) 

8 


The  dynamic  viscosity  of  the  gas  is  given  by  the  Sutherland  expression 

rj,  \  3/2 


P  =  PO 


/To  PS) 
T0J  VT  +  Pi 


(19) 


where  values  of  //() ,  P0  and  Pi  are  given  by  the  thermochemistry  module.  The  viscosities  for  mixtures  of 
gases  is  obtained  from  Wilke’s  [20]  expression 


pmix  — 


N 

E 

S=1 


f  S  Ps 


MW ,  ’ 


where 


N 


.s'=l 


$s  = 

MWS, 


1  + 


Ps 

Ps' 


1/2 


/  MW  s' 
V  MWS 


1/4 


1  + 


MWS 

MW* 


1-1/2 


Following  [13],  the  Prandtl  number  is  given  approximately  as 

207 


Pr  = 


397  -  15 


(20) 


(21) 


(22) 
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3.1.3  Data  Structures 

Data  for  each  slug  of  gas  is  stored  in  a  data  structure  (called  “slug.data”)  which  contains  an  array  of  “L_cell” 
structures  (one  element  for  each  Lagrangian  cell)  and  other  data  such  as  boundary  flags  and  time  step  in¬ 
formation.  Each  Lagrangian  cell  (L_cell)  structure  contains  the  location  of  the  cell  midpoint,  the  location 
of  the  interface  to  the  right,  the  mass  contained  by  the  cell,  a  cell  average  of  the  local  flow  state,  and  time 
derivatives  of  the  state  variables.  Full  details  may  be  obtained  from  the  header  file  “lld.h”. 


cell  i ncbcx 


ghost 


gWosfc  oc-fi  s 


5 

\vitcr-facx.  ind-nx  's  $ 

S 


Figure  8:  Data  storage  for  a  single  slug  of  gas.  The  indexing  for  the  cells  is  shown  above  the  array  while  the 
indexing  scheme  for  the  interfaces  is  shown  below. 


Figure  8  shows  the  indexing  arrangement  for  the  cells  within  each  slug  data  structure.  The  array  consists 
of  both  internal  cells  ( ixmin  <  ix  <  ixmax)  and  ghost  cells  at  each  end  ( ixmin  —  2 ,  ixmin  —  1 ,  ixmax  + 
1,  ixmax  +  2).  The  ghost  cells  lie  outside  the  physical  gas  slug  and  contain  data  used  in  the  application  of 
the  boundary  conditions. 


3.1.4  Internal-Interface  Pressures  and  Velocities 


The  pressures  and  velocities  used  in  equations  (3)  and  (4)  are  obtained  by  first  interpolating  the  flow  state 
(consisting  of  a  set  of  values  for  p,  u,  v.  e,  P,  a)  from  the  cell  centres  to  either  side  of  each  interface  at  the 
start  of  the  time  step  and  then  applying  a  Riemann  solver  to  estimate  the  flow  states  at  the  interfaces  during 
the  time  step. 

The  state  of  the  flow  either  side  of  each  interface  “L”  and  “R”  is  interpolated  (or  reconstructed)  from  the 
set  of  cell  averaged  states  by  assuming  a  linear  variation  of  the  variables  within  cells.  This  interpolation  is 
performed  separately  for  each  primary  variable.  For  example,  the  density  either  side  of  interface  {j  +  \)  is 
obtained  by  a  nonlinear  interpolation  (or  reconstruction)  using  the  expressions 


where 


PL  =  Pj  +  (xj+\  ~  xj)  MINMOD((A-)j,  (A+)j)  , 

PR  =  Pj+i  +  {xj+i  -  xj+i)  MINMOD((A-)i+i,  (A+)j+i)  , 


(A-), 


(A+ )j 


Pj  Pj-i 
xj  ~  xj- i 
Pj+ 1  —  Pj 


(23) 


(24) 
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represent  two  possible  estimates  of  the  slope  of  the  density  for  cell  j  and  Xj-\,  xj,  xJ+  \  are  the  midpoints  of 
the  cells  j  —  1,  j  and  j  +  1  respectively.  The  MINMOD  limiter  function  selects  the  slope  with  the  minimum 
magnitude  if  both  slopes  have  the  same  sign  and  returns  zero  otherwise  (see  e.g.  [21]). 

Interpolation  for  the  other  variables  is  done  similarly.  To  make  the  code  more  robust,  the  conditions 
Pl,Pr  >  Pm  in  and  e^,  eR  >  cmin  are  imposed  after  interpolation  but  before  the  application  of  the 
Riemann  solver.  Details  of  the  Riemann  solver  are  already  available  in  [22],  but  for  completeness  and 
because  the  solver  is  related  to  the  implementation  of  the  specified- velocity  boundary  condition,  a  complete 
description  is  included  here. 

The  Riemann  solver  used  here  is  a  2-stage  approximate  solver  in  which  the  first  stage  computes  the 
intermediate  pressure  and  velocity  assuming  isentropic  wave  interaction.  A  second  stage,  based  on  the 
strong-shock  relations,  may  be  invoked  to  improve  the  first-stage  estimate  if  the  pressure  jumps  across  either 
wave  are  sufficiently  large.  In  practice,  this  modification  has  been  required  only  in  extreme  conditions  [22] . 
If  stage  2  (strong  shock  modification)  is  not  invoked,  the  solver  is  much  like  Osher’s  approximate  Riemann 
solver  [23]. 


STAGE  1:  The  first  stage  of  the  Riemann  solver  assumes  that  a  spatially  constant  left  state  (subscript  L)  and 
right  state  (subscript  R)  interact  through  a  pair  of  finite-amplitude  (and  isentropic)  compression  or  rarefaction 
waves.  Perfect  gas  relations  ([24]  cited  in  [25])  are  used  to  obtain  the  intermediate  states  (L* .  R* )  in  the  gas 
after  the  passage  of  left-moving  and  right-moving  waves,  respectively.  The  expressions  implemented  in  the 
code  are 


Pl  =  P*=P*  =  PL 


{i-1){Ul~Ur) 
2ol(1  +  Z ) 


27/(7-!) 


(25) 


and 


where  the  Riemann  invariants  are 


=  UR  =  U 


UlZ  +  Ur 
1  +  Z 


(26) 


UL 


Ur 


UL  + 


2ol 

7-1 


UR  - 


2  or 
7-1 


and 


(27) 


and  the  intermediate  variable  Z  is  given  by 


Z= 


aL 


(7— 1)/(27) 


(28) 


In  the  exceptional  situation  of  {U  l  —  Ur)  <  0,  we  assume  that  a  (near)  vacuum  has  formed  at  the  cell 
interface  and  set  all  of  the  interface  quantities  to  minimum  values. 

STAGE  2:  If  the  pressure  jump  across  either  wave  is  large  (say,  a  factor  of  10),  then  the  guess  for  the 
intermediate  pressure  is  modified  using  the  strong  shock  relations. 

If  P*  >  10  Pl  and  P*  >  10  Pr  then  both  waves  are  taken  to  be  strong  shock  waves  and  the  intermediate 
pressure  and  velocity  can  be  determined  directly  as 


JD* 


7  +  1 


PL 


\[PR 


,\fPR  +  y/PL 


(' UL  ~  UR) 


(29) 


and 


*  _  s/PL  UL  +  y/PR  UR 

'  - ~r == - <  - 

y/PR  +  yfPL 


(30) 
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Figure  9:  Applying  boundary  conditions  via  ghost  cells:  (a)  reflecting  end  condition;  (b)  supersonic  outflow; 
(c)  gas-gas  interface  (or  direct  exchange  of  data). 


If  P*  is  greater  than  Pi  or  Pr  (but  not  both),  the  stage- 1  estimate  for  P*  can  be  improved  with  two 
Newton-Raphson  steps  of  the  form 


where 


and 


UL 


Ur 


D*  p*  p 

rn+ 1  rn  r  n 


dFn 

dP* 


-l 


Fn  =  U*L{P*n)-U*R{P* )  , 


TJr  _  2 at  | 

UL  7-1  l  Pl 


= 


= 


7~1 

27 


UL  - 

Ur  + 

UR  + 


IP 


Pl 


(7+1) ) 


\  1/2 


2gfl 

7-1 


7-1 

27 


P* 

Pr, 

(  2P*  \ 1/2 

\Pr( 7+1) ) 


P*  <  10  PL 
P*  >10  PL  . 

P*  <  10  PR 
P*  >  10  PR 


(31) 


(32) 


(33) 


(34) 


During  the  update,  we  ensure  that  P*  >  Pm  in  where  Pm  in  is  some  small  value.  After  updating  P* ,  the 
intermediate  velocity  is  evaluated  using  the  relevant  strong-shock  relation  from  (33)  or  (34). 

The  pressure  and  velocity  at  each  interface  may  now  be  substituted  back  into  equations  (1),  (3)  and  (4) 
to  give  the  motion  of  the  cell  interfaces  and  the  rate  of  change  of  momentum  and  energy  within  the  cells. 


3.1.5  Boundary  Conditions 

Before  interpolation,  the  inviscid  boundary  conditions  are  applied  by  setting  up  two  layers  of  ghost  cells 
along  each  of  the  boundaries.  This  is  shown  schematically  in  Fig.  9.  For  a  supersonic  inflow  boundary,  all 
of  the  ghost-cell  quantities  are  specified  as  fixed  while,  for  a  supersonic  outflow  boundary,  the  ghost-cell 
quantities  are  extrapolated  from  active  cells  just  inside  the  boundary.  Solid-wall  (i.e.  reflective)  boundary 
conditions  are  applied  by  setting  all  of  the  scalar  quantities  in  the  ghost  cells  equal  to  those  in  the  active  cells 
adjacent  to  the  boundary  but  setting  the  ghost-cell  velocities  to  the  negative  of  the  velocities  in  the  active 
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Figure  10:  Schematic  diagram  showing  the  pressure  forces  on  a  piston. 


cells.  Where  two  gas  slugs  interact,  data  between  the  two  end-cells  of  the  first  slug  and  the  corresponding 
ghost-cells  of  the  second  slug  are  exchanged  as  shown  in  Fig.  9. 

Where  the  gas  interacts  with  a  piston  (or  end  wall),  the  boundary-interface  velocity  u*  is  specified.  The 
interface  pressure  may  then  be  determined  from  the  isentropic  relations  which,  for  a  right-end  boundary,  give 


P * 


G Ul-u *) 


(7  - 1) 

2  71/2 


27/(7-!) 


(35) 


Similarly,  the  interface  pressure  at  a  left-end  boundary  is  given  as 


JD* 


(u* 


Ur) 


(7-1) 

27V2 


27/(7-!) 


(36) 


3.1.6  Piston  Dynamics 

Each  piston  is  assumed  to  have  fixed  mass  (rnp),  length  (Lp)  and  frontal  area  (Ap).  The  piston  state  is  given 
by  a  flag  indicating  whether  the  piston  is  constrained,  its  centroid  position  (xp)  and  its  velocity  (Vp).  The 
governing  differential  equations  are 

d  v 

dtxP  -  VP  , 

J,V-  =  £[*«.-'>)+*>]  .  (3^) 

where  Pr  and  Pp  are  the  pressures  on  the  “back”  and  “front”  piston  faces  respectively  and  Fj  is  the  total 
frictional  force.  Refer  to  Fig.  10  for  the  general  arrangement.  If  the  piston  is  initially  restrained,  a  specified 
value  of  back-face  pressure  (Pp)  must  be  exceeded  before  the  piston  is  released. 
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For  the  simulation  of  the  X2  facility,  the  frictional  force  is  assumed  to  be  due  to  the  “chevron”  seal  near 
the  front  face  of  the  piston.  The  maximum  magnitude  of  the  frictional  force  is 

\Ff\max  =  HfASealPF  >  (38) 

where  /jf  is  the  coefficient  of  friction  of  the  seal  material  on  the  tube  wall  (taken  to  be  0.2)  and  Aseai  is  the 
effective  frontal-area  of  the  seal.  The  actual  value  of  Ff  used  in  equation  (37)  is 

(  —sign(Vp)  \Ff\max  if  (\VP\  >  Vtoi  or \Ap(PB  -  Pf) \  >  \Ff\max) 

Ff  =  {  (39) 

l  ~Ap(Pb  -  Pf)  if  (\VP\  <  Vtoi  and  \Ap(PB  -  Pf) \  <  \Ff\max ) 

where  the  velocity  tolerance  is  Vtoi  =  10~6  m/s. 


3.1.7  Diaphragms 

Diaphragms  are  implemented  as  a  flag  for  the  status  of  the  diaphragm  (intact  or  burst)  and  a  burst  pressure. 
Note  that  the  burst  pressure  is  a  “dynamic”  burst  pressure  which  may  be  significantly  higher  than  the  burst 
pressure  obtained  in  hydrostatic  rupture  tests  [26].  The  effect  of  the  diaphragm  is  coded  directly  into  lld.cxx 
as  a  change  in  boundary  conditions  selected  by  the  diaphragm’s  status  flag.  For  example,  two  gas  slugs 
initially  separated  by  a  diaphragm  will  have  reflective  boundary  conditions  applied  at  the  diaphragm  station. 
On  rupture,  the  applied  boundary  conditions  will  be  changed  to  a  data-exchange  condition. 


3.1.8  Time  Stepping 

The  state  quantities  for  both  pistons  and  gas  slugs  are  advanced  from  time  level  n  to  time  level  n  +  1  with 
the  predictor-corrector  scheme 


A C/(1) 
t/U) 

A  U{2) 


At 


dU (") 


dt 

jjV)  +  Af7(1) 

dUM 


At 


dt 


U[1)  +  ^  [AU{2)  -  A 


(40) 


where  the  superscripts  (1)  and  (2)  indicate  intermediate  results  and  (^)  includes  the  rate  of  change  of 
interface  positions,  cell  momentum,  cell  energy,  piston  velocity  and  piston  position.  If  a  first-order  scheme 
is  desired,  only  the  first  stage  is  used  and  (/O'+F  =  [J(  l  K  Although  first-order  time-stepping  requires  fewer 
operations  than  second-order  time-stepping,  it  is  also  less  robust. 

To  maintain  stability,  the  time  step  is  restricted  to 


At  5:  Atanowec[  —  CFL  Atsignal  ,  (41) 

where  A tauowed  is  the  smallest  value  for  all  cells  (and  all  gas  slugs)  and  CFL  is  the  specified  Courant- 
Friedrichs-Lewy  number.  It  is  normally  restricted  to  CFL  <  0.5  in  the  simulations  discussed  later.  For  each 
cell,  the  inviscid  signal  time  is  approximated  as 


Af 


signal  — 


Ax 

ttl  +  a 


(42) 
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3.2  Axisymmetric  Flow  Code,  mbcns2 

For  simulating  two-dimensional  flow  fields,  we  use  the  code  mbcns2  [27]  is  an  extension  of  the  single¬ 
block  Navier-Stokes  integrator  CNS4U  [28]  to  structured  multiple-block  domains.  The  name  mbcns2  is  an 
acronym  for  Multiple-Block  Compressible  Navier-Stokes  solver,  second  version.  The  code  solves  the  com¬ 
pressible  Navier-Stokes  equations  via  a  cell-centred  time-dependent  finite-volume  formulation.  The  govern¬ 
ing  equations  are  expressed  in  integral  form  over  arbitrary  quadrilateral  cells,  with  the  time  rate  of  change  of 
conserved  quantities  in  each  cell  specified  as  a  summation  of  the  mass,  momentum  and  energy  flux  through 
the  cell  interfaces.  The  code  is  capable  of  considering  both  planar  and  axisymmetric  two-dimensional  ge¬ 
ometries  and  the  thermochemistry  module  can  handle  gases  in  chemical  equilibrium  or  nonequilibrium. 
When  simulating  gases  with  finite -rate  chemistry  and  radiation  energy  exchange,  these  physical  processes 
are  treated  with  an  operator-split  approach.  The  computational  core  of  mbcns2  is  written  in  a  combination 
of  C  and  C++,  with  the  option  for  user-defined  functions  such  as  boundary  conditions  provided  as  Lua  scripts. 
Preprocessing  ( i.e .  grid  generation)  and  postprocessing  is  handled  by  a  collection  of  Python  programs. 


3.2.1  Governing  Equations  for  the  Gas  Dynamics 


The  code  is  formulated  around  the  integral  form  of  the  Navier-Stokes  equations,  which  can  be  expressed  as 


d_ 

dt 


UdV 


'V 


■  h  dA  + 


(43) 


where  S  is  the  bounding  surface  and  h  is  the  outward-facing  unit  normal  of  the  control  surface.  For  axisym¬ 
metric  flow  as  considered  in  the  present  study,  V  is  the  volume  and  A  the  area  of  the  cell  boundary  per  unit 
radian  in  the  radial  direction. 


The  array  of  conserved  quantities  is  dependent  on  the  thermal  model  under  consideration,  and  for  the 
thermal  nonequilibrium  models  is 


U  = 


P 

pux 

puy 

pE 

PeVm 

pee 

Pfs  _ 


(44) 


Flere,  the  conserved  quantities  are  respectively  density,  .i-momcntum  per  volume,  y-momcntum  per  volume, 
total  energy  per  volume,  vibrational  energy  for  mode  m,  electronic-electron  energy  and  mass  density  of 
species  s.  Note  that  pee  includes  both  bound  and  free  electron  energy.  We  choose  to  solve  both  total  and  all 
individual  species  continuity  equations  to  add  rigour  to  our  solver:  the  redundant  information  gives  us  a  good 
idea  when  the  numerics  are  running  into  trouble.  Conversely,  when  only  solving  n  —  1  species  equations, 
it  is  easier  for  undetected  error  in  mass  fractions  to  accumulate.  Thus  for  11  species  air  with  6  vibrating 
molecules  and  the  inclusion  of  electrons,  for  example,  there  are  22  conserved  quantities. 


The  flux  vectors  are  divided  into  inviscid  and  viscous  contributions.  The  inviscid  component  in  thermal 
nonequilibrium  is 


pux 

pUy 

pu2x+p 

piixUy 

pUyUX 

pul  +  p 

pEux  +  pux 

i  + 

pEUy  +  pUy 

PeVm  U% 

Pevmuy 

pCeUx  Pe^x 

pCeUy  T  PeUy 

s 

Sr* 

*5 
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and  the  viscous  component  is 


Fv  = 


0 

0 

7~XX 

Txy 

Tyx 

Tyy 

T~xxUx  4“  PyxO-Ly  4“  qX 

i  + 

TxyTx  4“  Ty  y  Ty  4“  qy 

qx,vm 

%,Vm 

Qx,e 

Qy,e 

Jx,s 

Jy,s 

j  ■ 


(46) 


where  the  axisymmetric  viscous  stresses  are 


dux  ( dux  duy  Uy 

t~xx  —  2/r— - b  A  I  — - b  -x - 1 - 

ox  \  ox  dy  y 


duv 

Tyy  =  +  A 


dux  du,,  u. 


+ 


•y  _|_  ay 


dx  dy  y 


7 ~xy  Tyx  /i 


dux  dun 


(47) 


dy  dx 

and  where  the  secondary  viscosity  coefficient  A  is  expressed  in  terms  of  the  primary  coefficient  //  via  Stokes 

2 

'3> 


hypothesis,  A  =  —  1/r.  The  viscous  heat  fluxes  are 


qx  = 

qy  = 

qx,vm  = 
%,vm  = 
qx,e  = 

cly,e  = 


dT  ^  .  dTVn  ,  dTe 


ktr  Qx  +  E  kv 

s= mol. 


dx  "H"  dx 


+  E  Jx,shs  9 


s=all 


,  9T  \  -  dTVs  dTe  \  - 

ktr  Qy  +  E  +  fce  +  E 

s=mol.  s=all 


k  dTVm  ,  j  h 

/v?>™  ^  i  x.m1  H 


dx 


k  9TVm 

/vi  j ^  i 


dy 


y,m't'Vm  ’ 


rj  h 

e  ^  ,Jx,s'Les  > 

s=all 

A'  b  ^  ^  Jy,shes  • 

s= all 


dy 


(48) 


The  vector  of  source  terms  is  separated  into  geometric,  chemistry,  thermal  energy  exchange  and  radiation 
contributions  in  order  to  apply  the  operator-splitting  integration  approach,  Eq.  49. 


Q  —  Qgeom.  T  Qchem.  T  Q therm.  4~  Qrad. 
The  geometric  source  term  vector  for  axisymmetric  geometries  is 


(49) 


Q 


geom. 


o 

o 

(p  -  ree)  Axy/V 
0 
0 
0 


(50) 
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where  Axy  is  the  projected  area  of  the  cell  in  the  (x,y)-plane  and 


uv 

Tge  =  2  [i—  +  A 
y 


duy  .  uy\ 

dy  y  ) 


(51) 


For  planar  geometries  Qgeo m.  is  a  zero  vector.  See  the  original  ICASE  report  [28]  for  a  derivation  of  these 
terms. 

The  chemistry  source  term  vector  is 


Qchem. 


0 

0 

0 

0 


nvc 

Y  .  nEC 

A-is= ion.  J  Ls 


U> 


S 


(52) 


and  the  thermal  energy-exchange  source  term  vector  is 


Q  therm. 


0 

0 

0 

0 


0 


(53) 


The  radiation  source  term  vector  is 


rad.  — 


0 

0 

0 

V  •  Qrad 
0 

V  •  qrad 

0 


(54) 


where  any  purely  vibrational  component  of  radiative  heat  loss  (or  gain)  has  been  neglected.  The  transport, 
thermodynamic  and  chemical  kinetic  source  term  models  will  be  discussed  in  detail  in  Section  3.3. 


3.2.2  Discretised  Equations  and  Time-Stepping  Procedure 

The  conservation  equations  are  applied  to  straight-edged  quadrilateral  cells  for  which  the  boundary,  projected 
onto  the  (x,y)-plane,  consists  of  four  straight  lines.  These  lines  (or  cell  interfaces)  are  labelled  North,  East, 
South  and  West  and  the  integral  equation  is  approximated  as  the  algebraic  expression 


dU 

dt 


-y  'Yh  {Fi  ~  Fv)  ■  ndA  + Q  , 

V  NESW 


(55) 


where  U  and  Q  now  represent  cell-average  values.  An  operator- splitting  approach  as  advocated  by  Oran 
and  Boris  [29]  (see  Chapter  1 1  of  their  text)  is  applied  whereby  the  physical  mechanisms  are  applied  in  a 
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decoupled  fashion.  The  time  integration  of  the  ODE  system  shown  in  Eq.  55  is  then  approximated  by 


where, 


tt  ^  '  (-^i)  '  ft  dA  +  Qgeom.  T  Qvad.  > 

V  NESW 

{~Fv)-ndA, 

NESW 


chem. 


dt 


therm. 


Qchem.  ? 


Q therm.  • 


(56) 

(57) 

(58) 

(59) 

(60) 


Integration,  in  time,  of  the  discretised  equations  proceeds  in  a  loosely  coupled  fashion.  The  order  of  oper¬ 
ations  for  a  single  time-step  for  a  radiating  gas  in  thermochemical  nonequilibrium  is  shown  in  Figure  11. 
Some  of  the  chemical  kinetic  and  thermal  energy-exchange  ODE  systems  are  “stiff”  and  so  “subcycling” 
is  used  over  the  global  integration  time  step  via  smaller  steps  if  the  system  fails  to  solve.  The  number  of 
chemical  and  thermal  subcycles  are, 


Nc  =  At /Atc, 

Nt  =  At/Att. 

Currently  the  radiative  source  term  vector,  Qrad,  is  applied  closely  coupled  with  the  inviscid  fluxes.  This 
seems  to  be  adequate  for  the  work  done  thus  far,  but  may  need  to  be  revised  for  strongly  radiatively  coupled 
flows. 

The  advantage  of  the  operator-splitting  approach  is  that  the  optimal  integration  scheme  for  each  compo¬ 
nent  of  the  physics  can  be  implemented.  This  is  especially  useful  for  solving  large  chemical  kinetic  systems. 
The  resultant  set  of  ODE  systems  are  integrated  in  a  time  via  a  simple  predictor-corrector  scheme  for  the 
inviscid  and  viscous  increments,  one  of  a  selection  of  methods  (including  a  method  for  stiff  systems)  for 
the  chemistry  increment  (see  Section  3.3)  and  the  4th  order  Runge-Kutta-Fehlberg  method  for  the  thermal 
energy-exchange  increment. 

3.2.3  Multiple-Block  Grids,  Parallelisation  and  Boundary  Conditions 

As  shown  in  Figure  12,  the  data  arrays  for  each  block  are  dimensioned  such  that  there  is  a  buffer  region,  two 
cells  deep,  around  the  active  cells,  which  completely  defines  the  flow  domain  covered  by  the  block.  The 
buffer  region  contains  ghost  cells  which  are  used  to  hold  a  copy  of  the  flow  information  from  adjacent  blocks 
or  to  implement  the  boundary  conditions. 

For  a  boundary  common  to  two  blocks,  the  ghost  cells  in  the  buffer  region  of  each  block  overlap  the  active 
cells  of  the  adjacent  block.  The  only  interaction  that  occurs  between  blocks  is  the  exchange  of  boundary  data, 
prior  to  the  reconstruction  phase  of  each  time  step.  For  the  shared  memory  version  of  the  code,  the  exchange 
of  cell-average  data  along  the  block  boundaries  takes  place  as  a  direct  copy  from  the  active-cell  of  one  block 
to  the  ghost-cell  of  the  other  block.  Thus,  the  cells  along  the  common  boundary  of  each  block  must  match  in 
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1.  compute  gas  transport  due  to  inviscid  flux: 

(a)  apply  inviscid  boundary  conditions  or  exchange  data 
at  boundaries  for  each  block  as  appropriate 

(b)  reconstruct  the  flow  field  sate  on  both  sides  of  each  interface 

(c)  compute  the  inviscid  fluxes  Fi  •  n 

(d)  compute  the  radiative  source  term  —V  •  qrad  f°r  each  cell 

(e)  integrate  Eq.  57  over  the  timestep 

(f)  decode  the  conserved  quantities  via  an  equation-of-state  call 

(g)  repeat  for  corrector  update 

2.  compute  gas  transport  due  to  viscous  flux: 

(a)  apply  viscous  boundary  conditions  at  solid  walls 

(b)  compute  the  viscous  fluxes  as  Fv  •  n 

(c)  integrate  Eq.  58  over  the  timestep 

(d)  decode  the  conserved  quantities  via  an  equation-of-state  call 

(e)  repeat  for  corrector  update 

3.  compute  change  of  gas  state  due  to  chemical  reactions: 

(a)  compute  all  chemical  source  terms 

(b)  integrate  Eq.  59  over  the  timestep 

(c)  decode  the  conserved  quantities  via  an  equation-of-state  call 

(d)  redo  via  smaller  subcycles  if  failed  and  apply  call  to  equation-of-state  more  frequently 

4.  compute  change  of  gas  state  due  to  thermal  energy-exchange: 

(a)  compute  all  chemical  source  terms 

(b)  integrate  Eq.  60  over  the  timestep 

(c)  decode  the  conserved  quantities  via  an  equation-of-state  call 

(d)  redo  via  smaller  subcycles  if  failed  and  apply  call  to  equation-of-state  more  frequently 


Figure  11:  Sequence  of  operations  for  a  time-step  update  in  mbcns2. 
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Figure  12:  Active  and  ghost  cells  for  a  single  block  grid  for  mbcns2. 


both  number  and  position.  Some  logic  is  used  within  the  exchange  routines  to  set  the  appropriate  indexing 
direction  for  each  boundary.  The  information  on  the  connections  between  block  boundaries  is  stored  in  a 
(global)  connectivity  array.  For  each  boundary  on  each  block,  this  array  stores  the  identity  of  the  adjacent 
block  and  the  name  of  the  connecting  boundary  on  the  adjacent  block. 

Except  for  this  block  to  block  communication  (and  the  occasional  checking  of  time  step  magnitudes), 
the  rest  of  the  calculation  can  be  done  independently  for  all  blocks.  Thus,  the  algorithm  is  fairly  easy  to  im¬ 
plement  on  a  multiple-instruction,  multiple-data  (MIMD)  parallel  computer  and  we  have  a  single -program- 
multiple-data  (SPMD)  version  of  the  code  for  computationally-intensive  facility  calculations  as  shown  in 
Section  4.4.  When  running  such  simulations,  there  are  many  copies  of  the  program  running  independently 
on  separate  processors,  with  each  copy  of  the  program  handling  the  computation  for  a  single  block.  To  ex¬ 
change  block-boundary  data,  each  program  instance  must  communicate  with  the  other  programs  for  adjacent 
blocks.  The  communication  and  synchronisation  tasks  are  handled  via  a  standard  message  passing  library: 
MPI  [30], 

The  inviscid-component  of  other  boundary  conditions  is  also  applied  via  the  buffer  region  of  ghost  cells. 
Typically,  the  tunnel  surfaces  in  the  later  sections  are  assigned  fixed  temperature,  no-slip,  catalytic  (chemical 
equilibrium  at  wall  gas  state)  boundary  conditions.  Such  viscous  boundary  conditions  also  have  a  component 
implemented  using  data  at  the  cell  interfaces  that  lie  along  the  boundary  surface. 


3.2.4  Inviscid  Flux  Calculation 

The  flow-state  at  the  cell  interfaces  are  calculated  using  a  piecewise-parabolic  scheme  and  then  the  interface 
fluxes  can  then  be  calculated  via  an  appropriate  flux-calculator.  When  computing  the  inviscid  fluxes  at  each 
interface,  the  velocity  field  is  rotated  into  a  local  (n,p) -coordinate  system  with  unit  vectors 

h  =  nxi  +  ny  j  , 

P  =  Pxi  +  Pyj  ,  (61) 
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Figure  13:  Cell,  interface  and  vertex  indexing  for  mbcns2.  The  upper  half  of  the  figure  shows  the  primary 
cells  defining  the  finite- volumes  for  the  conservation  equations.  The  lower  part  of  the  figure  shows  the 
secondary  cells,  used  for  computing  spatial  derivatives. 
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normal  and  tangental  to  the  cell  interface  respectively.  We  have  chosen  the  tangential  direction  px  =  —ny 
and  py  =  nx.  The  normal  and  tangential  velocity  components 

Un  =  Tlx  Ux  T  Tly  Uy  , 

Up  —  Px  Ux  T  Py  Uy  ,  (62) 

are  then  used,  together  with  the  other  flow  properties  either  side  of  the  interface,  to  compute  the  fluxes 


r  f 

L  mass 

’  PUn 

F n— momentum 

punun  +  p 

Fp— momentum 

= 

punUp 

7? 

L  energy 

punE  +  pun 

F 

L  ±  species— s 

PUnfs 

in  the  local  reference  frame.  These  are  then  transformed  back  to  the  (x.  y)-planc  as 


Fmass 

Fmass 

F r— momentum 

F n— momentum^ x  F  F o—momentumPx 

F y— momentum 

= 

F n— momentum^ y  F  Fp—momenfurnPy 

F energy 

F energy 

F species— s 

F species— s 

For  the  simulation  of  shock  and  expansion  tubes,  the  shock  waves  can  be  extremely  strong  so  we  use 
the  default  adaptive  scheme  in  which  the  equilibrium  flux  method  (EFM)  [31]  is  applied  near  shocks  and  a 
modified  AUSMDV  calculator  [32]  is  applied  elsewhere. 


Reconstruction  The  primary  data  held  by  the  code  are  cell-average  data,  associated  with  cell  centres.  To 
get  the  fluxes  at  cell  interfaces,  a  variable-by- variable  reconstruction  is  made  of  the  flow  field.  This  is  done 
in  a  one-dimensional  fashion,  working  along  one-index  direction  at  a  time.  Left  and  Right  values  {wl  and 
wr  respectively)  of  a  flow  variable  at  a  cell  interface  arc  evaluated  as  the  corresponding  cell  average  value 
plus  a  limited  higher-order  interpolated  increment.  Given  an  array  of  cell-centres  [LI,  LO,  RO.  Ill]  with  an 
interface  located  between  LO  and  LO,  the  interpolated  values  are 


wl  =  wlo  +  [Al+  x  (2/?,£o  +  ^li)  +  A l -  x  Hro]  sl  , 
wr  =  wro  —  uro  [A#+  x  Ll o  +  A x  {2h.Ro  +  Lri)]  sr  , 


A  L-  = 

WLO  -  WL 1 

3  {hLO  +  h-Li) 

A  L+  = 

WRO  -  WLO 

3  {h-Ro  +  hL o) 

A  ^  ,  — 

wri  —  wRO 

^R+ 

3  (hm  +  h-Ri) 

olo  = 

hLo/Z 

h-Li  +  2/i£o  +  hLo 

«ito  = 

hRo/2 

hLo  +  Zh.Ro  +  hRi 

(65) 
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where  the  h  represents  the  width  of  a  cell  and  the  van  Albada  limiter  [33]  is  implemented  as 


Al Al+  +  |Al„Al+ 


A|_  +  A|+  +  e 
A/;-A^+  +  [Aj;_A 
A  r-  +  A  R+  +  e 


e 


1.0  x  10-12  . 


Finally,  minimum  and  maximum  limits  arc  applied  so  that  the  newly  interpolated  values  lie  within  the  range 
of  the  original  cell-centred  values.  Unlimited,  this  reconstruction  scheme  has  third-order  truncation  errors 
and,  with  the  limiter  active,  a  sine  function  is  reconstructed  with  an  effective  truncation  error  order  of  2.7. 

Typically,  reconstruction  is  done  for  density,  internal  energy,  velocity  components,  and  species  mass 
fractions.  Other  flow  quantities  that  are  needed  at  the  interfaces  for  the  inviscid  flux  calculation  are  then 
obtained  from  the  thermochemical  model. 


EFM  Calculation  The  equilibrium  flux  method  (EFM)  [31]  is  used  for  its  dissipative  nature  in  the  vicinity 
of  very  strong  compressions.  The  method  assumes  that  the  gas  is  in  equilibrium  and  the  molecular  velocities 
of  the  gas  either  side  of  the  interface  can  be  described  with  the  Boltzmann  distribution.  As  implemented  in 
Reference  [34],  the  flux  of  mass  from  the  left  state,  moving  to  the  right  is 


FmassL  =  W^  PL  UnL  +  Pl  \f2RTL  , 


(66) 


where 


er 


(67) 


Similarly,  the  flux  of  mass  from  the  right  state,  moving  to  the  left  is 


FmassR  =  WR  PR  UnR  +  DR  Pr  \J 2  RTr  , 


(68) 


where 


(69) 
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The  flux  vector  components  are  then 

Fmass  ^massL  T  F 'massR  ; 

Fn  —momentum  —  FmassL  R nL  F massR  R nR  PL  Pi?  > 

Fp— momentum  ~  FmassL  UpR  T  FmassR  UpR  , 


F 


energy 


=  (w+  PL  unL)  [eL  +  J-  +  * 

+  (^R  P-R  (  eR  +  “  F  -  ( U^R  +  UpR ) 

V  Pi?  Z  y 

+  (.Dj  \J2RTL  ppj  (u'L  +  UpL )  +  j- 

1  7  +  1 


+  (-Or  ^2RTr  pR^j  (u^r  +  UpR)  +  - 


FTr 


(70) 

v  /  y  z  *  z  7  —  ±  / 

Species  mass  fractions  are  just  transported  by  the  net  mass  flux  as  scalar  quantities.  Note  that  the  gas 
constants,  R  and  7,  are  not  really  constant;  they  are  density-weighted  averages  derived  from  the  local  values 
for  left  and  right  gas  states. 

AUSMDV  calculation  Most  of  the  flow  field  fluxes  are  computed  with  the  AUSMDV  [32]  because  of  its 
reasonably  low  dissipation.  The  calculation  procedure  starts  by  computing  the  weighting  parameters  for  the 
velocity  splitting 

2 Pl/pl 


«L  = 


OiR  = 


Pl/pl+Pr/pr 
2  Pr/ Pr 


Pl/pl+Pr/pr 

and  the  sound  speed  and  Mach  numbers  in  the  normal  direction  to  the  interface 


(71) 


&m  — 

rna x(aL,aR)  , 

Ml  = 

UnL 

^m 

Mr  = 

UnR 

(72) 

The  components  from  pressure  splitting  are  then 

PL, 


P+l  =  ^F(ML  +  iy  (2 -Ml) 


PL  duL 

R nL 
PR, 


,  otherwise  , 


pR  =  ^{Mr  -IY  (2  +  Mr) 


\Ml\  <  1.0  , 


\Mr\  <  1.0 


PR  duR 
UnR 


,  otherwise 


(73) 


where  duL  =  ^UnL+^UnL^  and  duR  =  (u'nR  The  components  from  the  velocity  splitting 

7  =  oil  ( ^ UnLA  +  - du^+duL  ,  \Ml\  <  1.0  , 


are 


„  4  ar 

=  duR  ,  otherwise  , 

( {UnR  CLrrr) 


uR  =  —  aR 


V  4  am 

=  duR  ,  otherwise  , 


+  duR  +  duR  ,  |  Mr  |  <  1.0  , 


(74) 
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These  components  are  then  combined  into  a  mass  flux 


(pu)  l  =U[  pL +  UR  pR 

and  a  pressure  flux 

Pi=Pl+  Pr 

and  a  normal-momentum  flux  (pu2)  1  as  a  blend  of  AUSMV  and  AUSMD  fluxes 

2 

(pu2) AUSMV  =  u~l  PL  UnL  +  U Pr  UnR  , 

(pu2) AUSMD  =  2  (( Pu)l{UuL  +  UnR )  ~  \  (pu)  1  \  (unR  —  , 

1  1 

(pu2)  1  =  (2  +  s)(pu2)AUSMV  +  (2  _  s)(pu2)ausmd  , 

with  the  switching  function,  s,  based  on  the  pressure  gradient 


1 

s  =  -  mm 
2 


\pr~pl\  ^ 
min  (pL,pR)J 


with  K  =  10. 

The  flux  vector  components  can  be  assembled  from  these  pieces  as 


(75) 

(76) 


(77) 


(78) 


7? 

- L  mass 
Fn— momentum 


(pu) 1  , 

(pu2)  1  +pi  , 
2  2 


(79) 


and  depending  on  which  way  the  wind  is  blowing,  the  remaining  flux  vector  components  are  assembled  from 
either  the  right  or  left  flow  properties.  For  (pu)  1  >  0, 


Fp— momentum 


Ff 


energy 


(pu)  l  UpL  , 
2 

(pu)  i  TTl  • 
2 


otherwise 


momentum 


energy 


(pu)  1  MpH  , 
2 


(80) 

(81) 


(82) 

(83) 


where  H  =  e  +  ^  +  1  (u^  +  u^)  is  the  total  enthalpy  of  the  gas.  Again,  species  mass  fractions  are  just 
transported  by  the  mass  flux  as  scalar  quantities. 

Finally,  an  entropy  fix  is  applied,  as  per  Section  3.5  in  Reference  [32].  This  first  determines  if  the  inter¬ 
face  includes  an  expansion  sonic  point 


Case  A:  unL  —  ol  <  0  and  unR  —  clr  >  0 
Case  B:  unR  +  ql  <  0  and  unR  +  or  >  0 
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and  increments  the  flux  if  only  a  single  expansion  wave  is  detected 


1  mass 


mass 


Aua  ( PR  ~  PL )  > 


F n— momentum 
Fp— momentum 


Ff 


energy 


^ ua  (pR^nR  PL^nL) 
^ ua  ( PRUpR  ~  PLUpL ) 
Aua  ( PrHr  -  PlHl) 


(84) 


where 


Aua  =  0.125((ttnR  —  an)  —  (unL  —  cll))  ,  for  A  and  not  B, 
=  0.125 ((unR  +  an)  -  (unL  +  aL))  ,  for  B  and  not  A, 
=  0  ,  otherwise. 


(85) 


Shock  Detector  The  switching  between  the  two  flux  calculators  is  governed  by  a  shock  (or  compression) 
detector.  This  is  simply  a  measure  of  the  relative  change  in  normal  velocity  at  interfaces.  Specifically,  we 
indicate  a  strong  compression  at  cell-interface  i  +  \  when 


(86) 


where  Tol  is  the  compression  tolerance  and  is  typically  set  at  -0.05.  This  measure  is  applied  to  all  interfaces 


in  a  block  and  then  a  second  pass  propagates  the  information  to  near-by  interfaces.  If  a  first  cell-interface  is 


identified  as  having  a  strong  compression,  the  EFM  flux  calculator  is  used  for  all  interfaces  attached  to  the 
cell  containing  that  first  cell-interface. 

3.2.5  Viscous  Fluxes 

The  viscous  flux  calculation  is  then  performed  based  on  the  the  updated  cell-centre  flow  state.  The  spatial 
derivates  required  in  the  viscous  stress  and  heat  flux  terms,  Eq.  47  and  48,  are  obtained  by  applying  the 
divergence  theorem  to  each  of  the  secondary  cells  surrounding  a  primary-cell  vertex  (as  shown  in  Figure  13, 
and  then  averaging  the  vertex  derivatives  to  obtain  a  cell  centre  value.  Secondary  cells  of  half  size  are  used 
along  the  boundaries  and  viscous  boundary  conditions  for  velocity  (e.g.  no  slip)  and  temperature  are  applied 
at  those  boundary  interfaces. 

3.3  Thermochemistry 

In  both  simulation  codes,  Lld2and  mbcns2,  the  governing  equations  are  closed  by  a  set  of  relations  be¬ 
tween  the  various  thermodynamic  properties  of  the  gas  mixture.  This  section  describes  a  gas  library  which 
provides  various  models  of  gas  behaviour  for  use  in  the  simulation  codes. 

There  are  a  number  of  gas  models  available  as  part  of  the  gas  library.  They  are  desribed  here: 

ideal  gas  mix 

The  ideal  gas  mix  is  used  to  model  one  or  more  components,  all  of  which  have  ideal  (perfectly  elastic 
collisions  and  calorically  perfect)  behaviour. 

equilibrium  gas  mix 

The  equilibrium  gas  mix  models  a  fixed-composition  gas  mix  which  is  assumed  to  be  in  thermal  and 
chemical  equilibrium  at  the  local  thermodynamic  conditions. 
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thermally  perfect  gas  mix 

The  thermally  perfect  gas  mix  models  a  gas  with  one  or  more  components,  all  of  which  have  perfect 
(collisional)  behaviour  but  each  have  all  internal  energy  modes  excited  to  an  equilibrium  described  by 
a  single  temperature 

multi-temperature  gas  mix 

This  model  represents  a  mixture  of  gases  with  perfect  collisional  behaviour  but  the  the  internal  en¬ 
ergy  modes  are  ascribed  different  Boltzmann  temperatures  in  order  to  model  the  effect  of  thermal 
nonequilibrium. 

Noble-Abel  gas  mix 

The  Nobel-Abel  gas  mix  model  is  used  in  interior  ballistics  work  where  high  pressures  mean  that  real 
gas  effects  cannot  be  neglected.  It  models  the  effect  of  the  finite  volume  occupied  by  gas  particles  when 
computing  the  equation  of  state.  Collisions  between  gas  particles  are  still  considered  to  be  perfectly 
elastic. 

van  der  Waals  gas  mix 

The  van  der  Waals  gas  mix  models  real  gas  effects  at  high  pressures/densities.  The  finite  volume 
occupied  by  the  gas  particles  and  the  non-ideal  particle  collisions  due  to  van  der  Waals’  forces  are 
considered  in  the  model. 

In  simulating  impulse  facilities,  we  typically  only  use  the  equilibrium  gas  and  thermally  perfect  gas 
models.  Therefore,  these  two  models  are  discussed  in  the  remainder  of  this  section.  The  thermally  perfect 
gas  mix  is  used  when  we  also  wish  to  simulate  the  effect  of  finite-rate  chemistry  between  the  gases.  The 
multi-temperature  gas  mix  has  also  been  used  recently  to  model  the  X2  facilities  [35],  however,  we  consider 
the  validation  of  this  model  a  work  in  progress. 

3.3.1  An  equilibrium  gas  mixture 

At  higher  temperatures,  the  gases  in  impulse  facilities  will  undergo  chemical  reactions.  For  example  in 
air,  by  2  000  K  oxygen  molecules  will  begin  to  dissociate,  and  by  4  000  K  nitrogen  molecules  will  begin  to 
dissociate.  When  the  chemical  reactions  are  rapid  compared  to  the  flow  transit  times,  we  use  a  model  of 
the  gas  which  assumes  chemical  equilibrium.  Along  with  the  assumption  of  chemical  equilibrium,  thermal 
equilibrium  is  also  assumed;  the  internal  energy  modes  of  the  gas  rapidly  equilibriate  with  the  translational 
temperature  when  compared  to  the  characteristic  flow  time.  This  assumption  is  most  appropriate  in  the 
high  pressure,  high  temperature  gas  that  is  driven  by  the  piston  against  an  unruptured  diaphgram:  the  high 
pressures  and  temperatures  give  rise  to  rapid  changes  in  chemical  composition. 

The  equilibrium  gas  mix  is  implemented  as  a  look-up  table  where  the  thermodynamic  properties  are 
interpolated  (or  extrapolated)  from  a  table  with  indexing  based  on  density  and  internal  energy.  The  look-up 
table  is  built  for  a  fixed  gas  composition  over  a  range  of  densities  and  energies  prior  to  the  gas  dynamics 
simulation  and  is  read  into  memory  at  the  start  of  the  simulation.  A  tool  provided  in  the  gas  library  builds 
the  look-up  table  by  running  the  CEA2  program  [36]  numerous  times. 

3.3.2  A  mixture  of  thermally  perfect  gases 

As  mentioned  above,  the  gas  model  for  a  mixture  of  thermally  perfect  gases  is  often  used  in  conjunction  with 
a  finite -rate  chemistry  simulation.  The  thermodynamic  relations  for  the  gas  mixture  are  presented  here.  The 
implementation  of  finite -rate  chemical  effects  is  discussed  in  Section  3.3.3. 


9-28 


RTO-EN-AVT-186 


CFD  Tools  for  Design  and  Simulation 
of  Transient  Flows  in  Hypersonic  Facilities 


A  single  thermally  perfect  gas  The  assumed  behaviour  of  a  thermally  perfect  gas  is  that  all  internal  energy 
modes  are  in  equilibrium  at  a  single  temperature.  For  atoms  this  means  that  the  Boltzmann  distributions  for 
translational  and  electronic  energy  are  governed  by  one  temperature  value.  Similarly  for  molecules,  the 
Boltzmann  distributions  for  translational,  rotational,  vibrational  and  electronic  energy  are  described  by  a 
single  temperature  value. 

To  model  a  thermally  perfect  gas  requires  a  knowledge  of  how  the  gas  stores  energy  as  a  function  of  tem¬ 
perature.  It  is  convenient  to  have  available  the  specific  heat  at  constant  pressure  as  a  function  of  temperature, 
CP(T).  From  this,  specific  enthalpy  of  the  gas  can  be  computed  as 


h 


CP{T)dT  +  h(Tref) 


ire/ 


and  entropy  is  given  as 


s  = 


Cp(T)  ,  nfrr 

dT  +  s(Tref). 


(87) 


(88) 


The  transport  properties,  viscosity  and  thermal  conductivity,  can  be  calculated  as  a  function  of  tempera¬ 
ture  for  a  single  component  of  the  gas  mix.  The  transport  properties  for  a  single  component  can  be  combined 
by  an  appropriate  mixing  rule  to  give  a  mixture  viscosity  and  thermal  conductivity. 

In  the  implementation,  a  thermally  perfect  gas  is  characterised  by  five  curve  fits  all  of  which  are  functions 
of  temperature: 

1.  specific  heat  at  constant  pressure,  CP(T), 

2.  enthalpy,  h(T), 

3.  entropy,  s(T), 

4.  viscosity, //(T),  and 

5.  thermal  conductivity,  k(T). 

The  form  of  these  curve  fits  follows  that  used  by  McBride  and  Gordon  [36].  The  curve  fits  for  thermodynamic 
properties  in  non-dimensional  form  are  as  follows: 

clqT  2  +  a\T  1  +  a2  +  cl^T  +  CI4T2  +  a^T3  +  clqT 4  (89) 

rri—2  1  m-ll  rp  I  ,  ^  T2  T3  T4  CL  7 

—a^T  +  a\T  logT  +  a2  +  <23— +  0:4  — — \- a§— — \~  o,q— — b  —  (90) 

2  6  4  hi 

rp—2  rp2  ^n4 

-a0-^ - a\T~l  +  a2  log  T  +  a3T  +  a4—  +  a5—  +  a6—  +  a8  (91) 

The  coefficients  for  these  curve  fits  are  available  for  a  large  number  of  gaseous  species  in  the  CEA  pro¬ 
gram  [36]  (and  associated  database  files).  Each  of  these  curve  fits  are  only  valid  over  a  limited  temperature 
range.  For  example,  the  thermodynamic  curve  fits  for  molecular  nitrogen  (N2)  are  comprised  of  three  seg¬ 
ments:  200.0-1000.0  K,  1000.0-6000.0  K  and  6000.0-20000.0  K.  Beyond  this  range  the  values  are  extrapo¬ 
lated  in  this  work.  The  extrapolations  are  based  on  a  crude  assumption  of  constant  Cp  outside  of  the  range. 


CP(T) 

R 

H(T) 
RT 
S{T ) 
R 
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Cp(Tiow) 

R 

Cp(Thjgh) 

R 

—  {H(Tiow)Tiow  -  C'p(TZow)(T/ol(;  -  T)} 

T  {H (Thigh)Thigh  +  Cp(Thigh)(T  —  Thigh )} 

S(Tlow)  -  Cp{Tlow)\og(^^j 

S (Thigh)  +  Cp(Thigh)  log  f  — - ^ 

V  1  high  J 

The  curve  fits  for  viscosity  and  thermal  conductivity  arc  also  in  the  same  form  as  that  used  by  the  CEA 
program  [36].  The  curves  are  as  follows. 

log  m(T)  =  a0  log  T  +  y  +  |§  +  a3 
log  k(T)  =  b0  log  T  +  ^  +  ||  +  63 

Mixing  rules  for  a  collection  of  thermally  perfect  gases  The  thermodynamic  state  for  a  mixture  of  ther¬ 
mally  perfect  gases  is  uniquely  defined  by  two  state  variables  and  the  mixture  composition.  The  internal 
energy  is  computed  as  a  mass  fraction  weighted  sum  of  individual  internal  energies, 

N  N 

e  =  =  fi  (bi  ~  RiT)  •  ^ 

i=l  i= 1 

Pressure  is  computed  from  Dalton’s  law  of  partial  pressures, 

N 

P  =  ^  PiRiT-  (93) 

2—1 

The  specific  gas  constant  for  the  mixture  is  defined  as 

N 

R  =  Y1  hRi-  (94) 

2—1 

The  calculation  of  Cp  is  based  on  a  mass  fraction  weighted  sum  of  component  specific  heats, 

N 

cP  =  Yl  f*cPi-  ^ 

i=  1 

The  specific  heat  at  constant  volume  is  then  computed  as 

Cv  =  Cp  —  R.  (96) 


Thus  the  extrapolations  are  as  follows: 

CP{T  <  Tlow) 

R 

Cp(T  >  Thigh ) 

R 

H(T  <  Tiow) 
RT 

H(T  >  Thigh ) 
RT 

S(T  <  Tlow) 

R 

S(T  >  Thigh) 
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The  ratio  of  specific  heats,  7,  is  given  by  its  definition, 


7  = 


Cp 

a, 


The  frozen  sound  speed  for  the  mixture,  a,  is  calculated  as 

a  =  y/'yRT. 


(97) 


(98) 


During  a  compressible  flow  simulation,  the  values  of  p  and  e  are  most  readily  available  from  the  con¬ 
served  quantities  that  arc  solved  for  during  each  time  increment.  This  leads  to  the  specific  problem  of  solving 
for  the  thermodynamic  state  of  the  gas  mixture  given  p,  e,  and  the  mixture  composition,  /  .  However,  the 
formulae  previously  presented  are  all  explicit  in  temperature.  We  solve  for  temperature  using  the  Newton 
iteration  technique  for  zero  solving, 

Tn+i  =  Tn-^\,  (99) 

J  0\ln) 

where  the  zero  function,  fo(T),  is  based  on  the  given  internal  energy,  e,  and  a  guess  for  internal  energy  based 
on  temperature, 

N 

MT )  —  &  guess  6  —  ^  ^  fi  (hi  Ri-^guess)  (100) 

i=l 

Using  the  fact  that  Cvi  =  ^,we  can  conveniently  find  the  derivative  function  for  the  Newton  technique  by 
computing  the  mixture  Cv, 

^  =  =  =  (io» 

i—  1  i— 1 

The  Newton  iteration  is  set  to  converge  when  the  accuracy  of  the  temperature  value  is  within  ±1.0  x  10-6  K. 
Personal  experience  has  shown  that  this  kind  of  error  tolerance  is  required  when  temperature  is  used  in  a 
finite -rate  chemistry  calculation  to  compute  rates  of  composition  change. 

The  calculation  of  mixture  transport  properties  is  not  as  straight  forward  as  the  thermodynamic  proper¬ 
ties.  A  mixing  rule  is  required  to  compute  the  mixture  viscosity  and  thermal  conductivity.  Wilke’s  mixing 
rule  [20]  has  been  implemented  in  the  work  presented  here.  Specifically,  the  mixing  rules  used  by  Gordon 
and  McBride  [37]  in  the  CEA  program  are  used  for  calculating  mixture  transport  properties  in  this  work; 
these  rules  are  a  variant  of  Wilke’s  original  formulation  [20]. 


N 


pmix  —  ^  ^ 


XiPi 


— ]  Xi  ±  y±,=i  Xj<f>ij 


and 


JL  r.u. 

j  _ 

Kmix  ~  ~~ 

i=l  xi  +  Li=!  xiVij 


(102) 


(103) 


where  x%  is  the  mole  fraction  of  species  i. 

The  interaction  potentials,  faj  and  iptj,  can  be  calculated  a  number  of  ways.  Again,  the  formulae  sug¬ 
gested  by  Gordon  and  McBride  [37]  have  been  used, 


^  =  I 


-(rori  ( 


2  M,- 


Mi  ±  M  j 


1/2 


(104) 
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and 

2.41(Mj  —  —  0.142Mj)" 

.  +  {Mi  +  Mj)2  . 

where  Mi  and  Mj  refer  to  the  molecular  weights  of  species  i  and  j  respectively. 


(105) 


Once  the  mixture  viscosity  and  thermal  conductivity  have  been  computed,  it  is  possible  to  compute  the 
mixture  Prandtl  number  from  its  definition 


Pr  = 


/iCp 

~k~ 


(106) 


3.3.3  Finite-rate  chemistry  implementation 

Rates  of  species  change  due  to  chemical  reaction  By  assuming  a  collection  of  simple  reversible  reactions, 
the  chemically  reacting  system  can  be  represented  as, 

N  N 

y;  ^Xi  ^  y;  PiXj,  (107) 

i=  1  2=1 

where  at  and  represent  the  stoichiometric  coefficients  for  the  reactants  and  products  respectively.  The 
case  of  an  irreversible  reaction  is  represented  by  setting  the  backward  rate  to  zero.  For  a  given  reaction  j,  the 
rate  of  concentration  change  of  species  i  is  given  as, 

=vi\kf  nw*  -  h  nro*  \  >  ^io8) 

3  l  i  i  ) 

where  vt  =  fit  —  ai.  By  summation  over  all  reactions,  Nr,  the  total  rate  of  concentration  change  is, 


d[Xj] 

dt 


iVy’ 

E 


3= 1 


d[Xj\ 

dt 


j 


(109) 


For  certain  integration  schemes  it  is  convenient  to  have  the  production  and  loss  rates  available  as  separate 
quantities.  In  this  case, 

Nr  Nr 

=  Qi  f'/  =  ^  '^apPi.j  ^  '  ^3yaj, ./  (HO) 

3= 1  i=1 

The  calculation  of  c oapPi  •  and  ujvai  .  depends  on  the  value  of  Vj  in  reach  reaction  j  as  shown  in  Table  1. 


Table  1 :  The  form  of  the  chemical  production  and  loss  terms  based  on  the  value  of  zy 


Vi  >  o 

o 

V 

S" 

Uappi  vikfY\i 

-VihWiiXiY* 

^ va,i  ^ikb 

v.kfU^Xir 

The  calculation  of  the  reaction  rate  coefficients,  kf  and  kb,  and  the  solution  methods  for  the  ordinary 
differential  equation  system  of  species  concentration  changes  are  discussed  in  the  subsequent  sections. 
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Reaction  rate  coefficients  The  reaction  rate  coefficients  for  a  reaction  can  be  determined  by  experiment 
(often  shock  tube  studies  are  used)  or  from  theory.  In  a  great  number  of  cases,  estimates  of  the  reaction  rate 
from  theory  can  vary  by  orders  of  magnitude  from  experimentally  determined  values.  For  this  reason,  fits  to 
experimental  values  are  most  commonly  used. 

For  the  single-temperature  gas  model  discussed  in  this  chapter,  the  forward  reaction  rate  coefficients  are 
calculated  using  the  generalised  Arrhenius  form, 

kf  =  ATnex  p(^)  (HD 

where  k  is  the  Boltzmann  constant  and  A,  n  and  Ea  are  constants  of  the  model. 

The  backward  rate  coefficient  can  also  be  calculated  using  a  modified  Arrhenius  form, 

kb  =  ATn  exp  (1 i2) 

or  it  can  be  calculated  by  first  calculating  the  equilibrium  constant  for  the  reaction, 

kb=^r.  (113) 

If  the  backward  rate  coefficient  is  calculated  from  the  equilibrium  constant,  then  a  method  of  calculation  of 
the  equilibrium  constant  is  required.  The  equilibrium  constant  for  a  specific  reaction  can  be  calculated  from 
curve  fits  or,  as  is  done  in  this  work,  using  the  principles  of  thermodynamics.  The  equilibrium  constant  based 
on  concentration  is  related  to  the  equilibrium  constant  based  on  pressure  by, 

<114) 

where  patm  is  atmospheric  pressure  in  Pascals,  7 7  is  the  universal  gas  constant,  v  =  u,  and 

Kp  =  exp  ■  (115) 

The  derivation  of  the  formula  for  Kp,  the  equilibrium  constant  based  on  partial  pressures,  can  be  found  in  any 
introductory  text  on  classical  thermodynamics  which  covers  chemical  equilibrium.  The  differential  Gibbs 
function  for  the  reaction,  AG,  is  calculated  using 


Ns 

A  G  =  s£JViGi  (H6) 

i 

where  each  G,  is  computed  from  the  definition  of  Gibbs  free  energy, 

Gi(T)  =  Hi(T)  —  T  x  Si(T)  (117) 

and  Gr  is  in  units  of  J/mol.  Hi  and  St  can  be  computed  in  the  appropriate  units  by  using  the  CEA  polynomials 
and  multiplying  by  77T  and  77  respectively. 

Some  caution  should  be  exercised  in  the  selection  and  use  of  reaction  rates  for  a  specific  flow  problem. 
In  many  cases,  a  set  of  reaction  rates  may  only  be  “tuned”  for  a  specific  problem  domain.  This  problem  of 
“tuned”  sets  of  reaction  rates  and  an  explanation  for  why  it  arises  is  described  by  Oran  and  Boris  (p.  38  of 
Ref.  [38]): 
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A  problem  that  often  arises  in  chemical  reactions  is  that  there  are  fundamental  inconsistencies 
in  a  measured  reaction  rate.  For  example,  there  may  be  experimental  measurements  of  both  the 
forward  and  reverse  rate  constants,  kf  and  kr.  Nonetheless,  when  either  is  combined  with  the 
equilibrium  coefficient  for  that  reaction,  the  other  is  not  produced.  This  appears  to  represent  a 
violation  of  equilibrium  thermodynamics.  The  explanation  is  usually  that  kf  and  kr  have  been 
measured  at  rather  different  temperatures  or  pressures,  and  so  there  are  inconsistencies  when 
they  are  extrapolated  outside  the  regime  of  validity  of  the  experiments. 

Solving  the  chemical  kinetic  ordinary  differential  equation  The  system  represented  in  Equation  109 
is  a  system  of  ordinary  differential  equations  (ODEs)  which  can  be  solved  by  an  appropriate  method.  For 
certain  chemical  systems,  the  governing  ODEs  form  a  stiff  system  due  to  rates  of  change  varying  by  orders 
of  magnitude  for  certain  species.  For  these  systems,  special  methods  for  stiff  ODEs  are  required.  In  this 
work,  four  methods  for  the  numerical  solution  of  the  ODE  system  have  been  implemented. 

1 .  Euler  method 

2.  modified  Euler  method 

3.  alpha-QSS  method,  and 

4.  Runge-Kutta-Fehlberg  method 

The  Euler  method  and  modified  Euler  method  are  standard  techniques  for  solving  ODEs  and  the  details 
can  be  found  in  any  text  dealing  with  numerical  methods  and  numerical  analysis.  The  fourth-order  Runge- 
Kutta  method  uses  a  fifth-order  error  estimate  as  a  means  for  controlling  the  timestep  used  for  integration  as 
proposed  by  Fehlberg  [39],  This  is  particularly  efficient  for  non-stiff  systems. 

alpha-QSS  method  The  alpha-QSS  (quasi-steady-state)  method  was  proposed  in  Mott’s  thesis  [40].  It 
is  an  ODE  solver  aimed  specifically  at  the  problem  of  stiffness  in  chemical  systems.  This  ODE  solver 
makes  use  of  the  forward  and  backwards  rates  of  concentration  change  as  calculated  by  Equation  110.  This 
is  a  predictor-corrector  type  scheme  in  which  the  corrector  is  iterated  upon  until  a  desired  convergence  is 
achieved.  The  predictor  and  corrector  are. 


(118) 


[Xi]n+1 


(119) 


In  the  above  equations, 


(4?  +  £?) 


(120) 


and 


Qi  =  dig?  +  (1  -  cti)  Q°. 


(121) 


The  key  to  the  scheme  is  calculating  a  correctly.  This  a  parameter  controls  how  the  update  works  on  a  given 
species  integration.  Note  that  a  is  defined  as 


(122) 
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Using  Pade’s  approximation, 


360  +  120x+  12x2 
360^-^403T+T2x^^1^x3~+~^ 


(123) 


it  is  possible  to  write  a  form  of  the  expression  for  a  which  is  more  amenable  to  computation  as  the  expensive 
exponential  function  evaluation  is  avoided.  The  approximation  for  a  becomes, 


.  ,  180r3  +  60r2  +  11?’  +  1 

^  ~  360r3  +  60r2  +  12r  +  1 


(124) 


where  r  =  l/(LAt). 


Coupling  chemistry  effects  to  the  flow  solver  Some  details  about  the  coupling  of  the  chemistry  effects  to 
the  gas  dynamics  simulation  are  provided  here.  In  an  unsteady,  time-accurate  flow  simulation,  the  allowable 
timestep  is  constrained  by  the  Courant-Friedrichs-Lewy  (CFL)  criterion.  In  a  viscous  compressible  flow,  the 
CFL  criterion  allows  one  to  select  an  appropriate  timestep  and  limit  the  propagation  of  flow  information  to 
distances  less  than  one  cell-width.  The  speed  at  which  flow  information  propagates  is  a  function  of  inviscid 
wave  speeds  and  viscous  effects. 

When  the  effects  of  finite -rate  chemistry  are  ‘split’  from  the  flow  simulation,  the  chemical  update  is 
solved  in  a  separate  step  in  which  the  flow  is  held  frozen.  (In  fact,  in  true  timestep-splitting,  all  other 
contributing  physics  is  frozen  during  the  chemistry  update.)  Thus  the  chemistry  problem  is  to  find  the 
updated  species  composition  at  the  end  of  the  flow  timestep. 

It  may  be,  and  is  quite  likely,  that  the  flow  timestep  is  not  an  appropriate  timestep  to  solve  the  chemical 
kinetic  ODE  problem.  When  the  timestep  for  the  chemistry  problem  is  smaller  than  the  flow  timestep, 
the  chemistry  problem  is  subcycled  a  number  of  times  until  the  total  elapsed  time  equals  that  of  the  flow 
timestep.  It  is  common  to  have  simulations  where  the  chemistry  timestep  is  100-1000  times  smaller  than  the 
flow  timestep,  that  is,  100-1000  subcycles  are  required  to  solve  the  chemistry  problem.  When  the  timstep  for 
the  chemistry  problem  is  larger  than  the  flow  timestep,  it  is  simply  set  to  the  value  of  the  flow  timestep. 

During  the  simulation  process,  the  chemistry  timestep  is  tracked  for  each  finite- volume  cell  in  the  simu¬ 
lation.  Although  the  flow  ‘moves  on’  in  subsequent  timesteps,  if  the  change  of  flow  conditions  is  not  large, 
then  the  previous  chemistry  timestep  will  be  a  good  estimate  to  begin  the  new  chemistry  problem  in  the  sub¬ 
sequent  timestep.  An  exceptional  case  is  when  a  shock  passes  through  the  cell:  the  change  of  flow  conditions 
does  become  large.  In  this  instance,  the  old  chemistry  step  is  disregarded  and  a  new  step  is  selected.  The 
selection  procedure  for  a  new  step  is  discussed  in  the  next  paragraph.  When  using  either  the  Runge-Kutta- 
Fehlberg  or  the  alpha-QSS  methods,  an  estimate  of  the  new  chemistry  timestep  is  provided  as  part  of  the 
ODE  update  routine. 

So,  during  a  simulation,  the  old  chemistry  step  at  one  iteration  is  used  to  begin  the  new  chemistry  problem 
in  the  next  iteration.  What  is  needed  is  a  means  to  select  the  chemistry  step  on  the  initial  iteration,  or 
whenever  the  old  suggestion  is  not  reasonable  (as  in  the  case  of  a  shock  passing  through  the  cell).  In  this 
work,  the  initial  step  for  the  chemistry  problem  is  selected  based  on  the  suggestion  by  Young  and  Boris  [41], 


dtchem  =  U  min 


axm\ 

\[*m) 


(125) 


where  e\  is  taken  as  1.0  x  HE  in  this  work,  and  the  expression  is  evaulated  at  the  initial  values  for  the 
chemistry  subproblem.  Young  and  Boris  [41]  suggest  that  e\  be  scaled  from  the  convergence  criteria.  We 
have  found  that  the  fixed  value  is  adequate  for  the  problems  of  interest  to  our  research  group. 
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4  APPLICATIONS 


In  Section  4.1,  the  applications  start  with  the  simplest  model  of  the  smaller  machine,  the  Drummond  tun¬ 
nel.  The  Lld2  simulation  gives  quick  results  for  bulk  gas  conditions  in  the  shock  tube  but  doesn’t  include 
the  details  of  the  reflected-shock  interaction  with  the  shock-tube  boundary-layer.  To  capture  these  details, 
and  the  expansion  process  through  a  high  Mach  number  nozzle,  axisymmetric  flow  simulations  are  used 
(Section  4.2).  Even  though  the  Drummond  tunnels  is  a  relatively  low  enthalpy  machine,  the  axisymmetric 
simulations  of  the  full  facility  are  computationally  demanding. 

For  the  modelling  of  the  X2  expansion  tube,  we  again  start  with  the  L 1  d2  modelling  of  the  whole  facility 
to  look  at  developing  operating  conditions  with  a  new,  light  weight  piston  (Section  4.3).  Since  this  mode  of 
analysis  in  exploratory,  we  need  to  do  many  simulations  and  concentrate  on  the  large-scale  facility  behaviours 
of  piston  motion  and  shock  speed.  Once  the  operational  conditions  are  defined,  we  need  to  increase  the  level 
of  detail  of  the  simulation  model  in  order  to  get  complete  estimates  of  the  test-section  flow  conditions.  The 
interaction  of  the  final  expansion  and  chemistry,  together  with  substantial  boundary  layers  growing  in  the 
acceleration  tube,  demands  a  high  resolution  numerical  simulation.  The  computational  cost  for  doing  this 
high-resolution  axisymmetric  simulation  for  the  entire  X2  machine  is  very  high  so  we  make  use  of  both 
simulation  codes.  The  Lld2  code  is  used  to  model  the  whole  facility  and  transient  flow  data,  just  upstream 
of  the  secondary  diaphragm,  is  fed  to  axisymmetric  simulation  of  the  unsteady  expansion  process  through 
the  acceleration  tube  and  into  the  dump  tank.  An  example  of  this  hybrid  simulation  approach  is  discussed  in 
Section  4.4. 

4.1  One-Dimensional  Modelling  of  the  Drummond  Tunnel 

Using  a  nitrogen  driver  gas  with  nitrogen  test  gas  gives  a  limited  but  simple-to-model  operating  conditon  for 
the  Drummond  tunnel.  The  Lld2  input  listing,  shown  below,  defines  the  simulation  by  specifying: 

•  the  gas  model, 

•  the  internal  surface  description  of  the  tunnel  as  a  set  of  break-points  in  diameter,  and 

•  the  gas  slugs,  which  are  contained  by  the  tube  end  walls,  a  contact  surface  (where  the  main  diaphragm 
would  have  been)  and  a  diaphragm  at  the  nozzle  throat. 

The  sudden  contractions  or  expansions  in  tube  diameter  at  the  diaphragm  and  nozzle  throat  are  also  marked 
as  regions  of  head  loss.  Figure  14  shows  the  representation  of  the  internal  surfaces  of  the  machine.  Note  that 
all  area  transitions  occur  over  finite  axial  distances.  This  is  a  limitation  of  the  quasi-one-dimensional  model 
implemented  by  Lld2,  but  also  happens  to  be  close  to  reality  for  the  Drummond  tunnel.  Instrumentation 
locations  in  the  actual  machine  are  identified  as  “history  locations”  in  the  listing  and  data  are  written  out  for 
those  locations  at  a  high  frequency  (every  2  qs).  Data  for  the  entire  flow  path  is  written  out  every  30  //s. 


#  "/ cf cfd2/app/Lld/examples/Drummond/dn2 .py 
gdata. title  =  \ 

"Drummond  tunnel  M4  nozzle  P4  =  3.25MPa  N2,  PI  =  30kPa  N2,  06-Jul-05 
gdata . case_id  =  27 
gdata . gas_name  =  "LUT" 

#  Define  the  tube  walls, 
gdata. n  =  4000 

add_break_point (-3 . 785,  0.0585)  #  upstream-end  of  the  driver  tube 

add_break_point (-3 . 035,  0.0585) 

add_break_point (-3 . 015,  0.0620)  #  steel-diaphragm  station 

#  Note  that  there  is  no  steel  diaphragm  in  this  simulation. 

add_break_point (  0.000,  0.0620)  #  downstrem  end  of  shock  tube 

add_break_point (  0.043,  0.0620)  #  start  of  contraction  to  throat 
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add_break_point (  0.080,  0.0220)  #  start  of  throat 

add_break_point (  0.100,  0.0220,  1)  #  start  of  nozzle  (conical  shape) 

add_break_point (  0.2653,  0.0700,  1)  #  start  of  parallel  test  section 

add_break_point (  0.30,  0.0700) 

# 

add_loss_region (-3 . 050,  -3.000,  0.5)  #  at  steel-diaphragm  station 
add_loss_region (  0.050,  0.120,  0.5)  #  at  nozzle  throat 
# 

gdata . T_nominal  =  296.0 


#  Create  the  gas-path. 

left_end  =  VelocityEnd (x0=-3 . 785,  v=0.0) 

driver_gas  =  GasSlug (p=3 . 25e6,  u=0.0,  T=296.0,  nn=150, 

viscous_ef fects=l ,  hcells=l,  label="driver  gas") 
cs  =  Gaslnterface (x0=-3 . 015) 

test_gas  =  GasSlug (p=30 . 0e3,  u=0.0,  T=296.0,  nn=300, 

viscous_ef fects=l,  hcells=l,  label="test  gas") 
diaph  =  Diaphragm (x0=0 . 1 0 ,  p_burst=150 . 0e3 ) 

dump_tank_gas  =  GasSlug (p=4 . 0e3,  u=0.0,  T=296.0,  nn=6,  to_end_L=l,  cluster_strength=l . 1 , 

viscous_ef f ects=l ,  hcells=l,  label="dump-tank  gas") 
right_end  =  FreeEnd (x0=0 . 3) 

assemble_gas_path (left_end,  driver_gas,  cs,  test_gas,  diaph,  dump_tank_gas,  right_end) 


#  Set  some  time-stepping  parameters 
gdata . dt_init  =  0.5e-6 
gdata. cfl  =  0.25 
gdata .max_time  =  8 . 0e-3 
gdata .max_step  =  100000 
add_dt_plot (0 . 0,  30.0e-6, 
add_history_loc (-0 . 2 95)  # 
add_history_loc (-0 . 078)  # 
add_history_loc (  0.000)  # 
add_history_loc (  0.090)  # 
add_history_loc (  0.265)  # 


2 . Oe-6) 

0,  heat  flux  gauge 
1,  pressure  transducer 
joint  at  nozzle  block 
mid-point  of  nozzle  throat 
nozzle  exit  plane 


2, 

3, 

4, 


Drummond  Tunnel  Profile  Drummond  Tunnel  Nozzle  Profile 


Figure  14:  Drummond  tunnel  cross-section  as  modelled  in  Lld2. 


From  the  data  recorded  throughout  the  simulation,  the  pressure  values  have  been  contoured  on  a  logarith¬ 
mic  scale  to  form  the  space-time  diagram  shown  in  Figure  15.  This  clearly  shows  the  shock  and  expansion 
wave  motion  within  the  machine,  including  the  shock  reflection  at  t  =  4  ms  and  over-tailored  interaction  of 
the  reflected  shock  with  the  test-gas  driver-gas  interface  at  t  =  5  ms,  both  occuring  at  the  downstream-end 
of  the  shock  tube. 

Comparison  with  experimental  data  is  shown  in  Figure  16.  There  is  good  agreement,  however,  the 
simulation  shows  a  number  of  sudden  jumps  in  pressure,  whilst  the  experimental  traces  show  more  gradual 
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ORGANIZATION 


dn.2_log_p.gen  logP.Pa  Space-time-plot 

xl  x2  dx  — 4.00e+00  5.00e-01  1.00e+00  yl  y2  dy  0.00e+00  8.00e-03  1.00e-03 
vl  v2  dv  2.74e+00  6.39e+00  2.43e-01 

O 

O 


Figure  15:  Computed  wave  diagram  for  Drummond  tunnel  operated  with  nitrogen  driving  nitrogen,  over¬ 
tailored  operation.  This  figure  is  scanned  from  PJ’s  workbook  and  the  annotation  relates  to  the  set  up  of  the 
axisymmetric  simulation  discussed  in  the  following  section. 
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changes  at  later  times.  This  is  a  limitation  of  the  quasi-one-dimensional  nature  of  the  simulation  and  the 
following  section  will  show  the  far  more  interesting  multidimensional  wave  interactions  that  occur  after 
shock  reflection. 


a)  n  ozzte  rese  rvotr  p  re  ssu  re  flow  lota!  temperature 


Figure  16:  Nozzle-supply  pressure,  Pitot  pressure  and  total  temperature  for  Drummond  tunnel,  nitrogen 
driving  nitrogen,  over-tailored  operation. 


Comparison  with  experimental  data  was  also  made  in  Reference  [42].  In  this  case,  the  pressure  of  the 
helium  driver  gas  at  diaphragm  rupture  was  3.25  MPa  and  the  filling  pressure  of  the  nitrogen  in  the  shock 
tube  was  set  at  16.5  kPa  to  achieve  tailored  operation.  Nozzle-supply  pressure  and  Pitot  pressure  histories 
are  shown  in  Figure  17.  The  shock  speed,  as  estimated  by  Lld2,  was  1.49  km/s  and  the  nozzle  supply 
conditions  were  ps  =  2.14  MPa  and  Ts  =  1920  K.  The  actual  (measured)  nozzle-supply  pressure  was 
consistently  around  2.0  MPa  for  a  number  of  shots.  These  conditions  are  very  similar  to  those  produced  with 
an  earlier  mode  of  operation  using  a  double-diaphragm  to  control  the  initial  release  of  the  driver  gas. 

The  simulation  has  accurately  captured  the  incident  shock  features,  both  shock  speed  and  pressure  jump. 
It  has  also  accurately  captured  the  long  term  decay  of  the  nozzle-supply  pressure  due  to  the  relatively  short 
driver  tube.  Where  the  simulation  deviates  noticeably  is  on  shock  reflection  where  the  gas  dynamics  is  quite 
complex.  Here,  Lld2’s  flow  approximation  using  only  normal  waves  is  lacking  in  the  detail  that  can  be 
captured  with  an  axisymmetric  analysis  using  mbcns2. 

4.2  Axisymmetric  Simulation  of  the  Drummond  Tunnel 

Starting  Flow  in  the  Mach  7  Nozzle  The  next  level  of  detail  is  captured  by  using  mbcns2  to  do  an 
axisymmetric-flow  simulation  of  the  last  0.5  metre  section  of  the  shock  tube  and  the  attached  supersonic 
nozzle.  Figure  18  shows  the  evolution  of  the  shock  reflection  and  nozzle-starting  process  for  the  Mach  7 
contoured  nozzle  (designed  by  Craddock  [2]). 

In  this  simulation,  the  gas  in  the  shock  tube  is  ideal  nitrogen,  initially  at  pi  =  30kPa  and  Tj  =  296  K. 
With  an  incident  shock  speed  of  782  m/s  from  the  Lld2  simulation,  the  post-shock  conditions  applied  at  the 
inflow  plane  are  p2  =  176  kPa,  T2  =  564  K  and  ug  =  535  m/s.  These  conditions  are  applied  to  the  core  flow 
of  the  inlet  and  a  boundary  layer  velocity  profile  is  superimposed.  The  simulation  starts  with  the  incident 
shock  at  the  inflow  plane,  and  the  boundary  layer  profile  grows  in  time. 

The  first  frame  in  Figure  18  occurs  shortly  after  the  reflection  of  the  shock  off  the  tube  endwall  and 
the  frames  continue  until  the  flow  through  the  nozzle  has  mostly  settled,  although  small  distrubances  are 
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Figure  1 7 :  Wall-pressure  history  in  the  shock-reflection  region  and  Pitot-pressure  history  at  the  nozzle  exit 
plane  for  single-diaphragm  operation  (shot  number  240397a3). 


still  working  their  way  through  the  nozzle.  There  is  much  detail  in  the  shock-reflection  region  but  the  main 
interest  of  this  simulation  is  the  complex  starting  process  in  the  expanding  part  of  the  nozzle.  The  classic, 
nearly  one-dimensional  starting  structure  [43]  seen  in  the  first  frame  evolves  into  an  extended  structure  with 
large  separation  zones  supporting  the  oblique  upstream-facing  shocks.  These  shocks  cause  recompression 
that  focuses  towards  the  nozzle  centreline.  Because  the  separation  line  takes  a  relatively  long  time  to  be 
swept  out  of  the  nozzle,  the  disturbance  near  the  centreline  persists  long  after  the  one-dimensional  theory 
indicates  that  the  flow  should  have  settled. 

Reflected-Shock  Interaction  with  the  Boundary  Layer  The  next  step  up  in  computational  effort  is  an 
axisymmetric  simulation  of  the  whole  machine.  Figure  19  shows  the  evolution  of  the  flow  in  the  shock- 
reflection  region  for  the  over-tailored  case  of  nitrogen  driving  nitrogen  [3].  There  is  much  of  interest  to  be 
observed  and  there  are  precious  few  one-dimensional  waves  involved.  The  idea  that  the  test  gas  is  brought 
to  rest  by  a  reflected  shock  wave  is  seen  to  be  far  from  the  more  detailed  reality. 

The  main  modelling  complexity  in  this  simulation  is  for  the  relatively  slow  opening  process  of  the  main 
diaphragm.  Here  the  metal  diaphragm  is  approximated  as  an  opening  iris,  with  the  barrier  between  driver 
gas  and  test  gas  being  removed  from  the  computational  domain  over  a  number  of  time  steps.  The  rest  of  the 
flow  domain  is  bounded  by  solid  walls,  at  296  K,  with  no-slip  boundary  conditions. 
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t  =  0.470  ms 


t  =  1.070  ms 


Figure  18:  The  development  of  flow  in  the  Mach  7  nozzle  of  Drummond  tunnel.  The  coloured  contours 
indicate  temperatures,  with  red  indicating  up  to  2600  K  and  blue  as  low  as  100  K. 
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Figure  19:  Sequence  of  numerical  schlieren  images  (top)  and  driver  gas  mass  fractions  ( bottom )  showing  the 
over-tailored  interaction  of  the  reflected  shock  with  the  boundary  layer. 
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4.3  Commissioning  a  New  Lightweight  Piston  for  X2 

4.3.1  The  Free  Piston  Driver 

The  free  piston  driver  is  currently  our  preferred  technique  to  achieve  high  performance  drivers  for  our  impulse 
facilities.  Since  this  type  of  driver  compressively  heats  the  driver  gas,  it  is  capable  of  achieving  both  the  high 
pressures  and  the  high  sound  speeds  required  to  generate  strong  shock  waves.  However,  the  free  piston  driver 
presents  three  main  challenges: 

1.  Tuning  its  operation  to  prevent  damage  to  the  piston,  to  the  buffer  and  to  the  compression  tube. 

2.  Tuning  its  operation  to  achieve  sufficient  constancy  of  driver  pressure  for  a  sufficient  duration. 

3.  The  monetary  cost  of  the  driver  assembly  is  very  high  since  the  structural  requirements  are  very  high. 


In  order  to  achieve  a  high  temperature  in  the  driver  gas,  large  compression  ratios  are  typically  used, 
therefore  the  volume  of  driver  gas  at  diaphragm  rupture  is  relatively  small.  If  the  piston  is  moving  with 
relatively  low  velocity  at  this  point,  the  driver  gas  slug  has  approximately  constant  volume.  The  unsteady 
expansion  will  therefore  lead  to  a  rapid  pressure  drop  in  the  driver  gas.  The  effect  of  this  pressure  drop  is 
then  transmitted  downstream  as  a  reflected  u  +  a  characteristic,  potentially  interfering  with  downstream  flow 
processes  before  or  during  the  test  time  [44]. 

We  have  recently  been  trying  to  produce  high  Mach  number,  high  static  pressure  flow  conditions  in  the 
X2  expansion  tube  facility,  however,  initial  attempts  did  not  achieve  expected  results.  The  existing  35  kg 
piston  is  relatively  heavy  for  the  length  of  compression  tube  and  therefore  is  operated  at  slow  speeds;  the 
result  being  that  the  driver  gas  maintains  its  pressure  for  a  relatively  short  duration.  For  the  high  speed  flow 
conditions  for  which  X2  is  typically  used  (such  as  planetary  entry  between  6  and  10  km/s),  critical  flow 
processes  occur  in  the  test  section  before  the  reflected  u  +  a  characteristic  from  the  driver  reaches  the  test 
section,  and  target  flow  conditions  are  therefore  achieved. 

However,  considering  the  slower  shocks  generated  through  the  dense  test  gas  for  high  total  pressure 
conditions,  early  pressure  loss  in  the  driver  manifests  itself  in  shock  speeds  which  rapidly  slow  down  before 
the  critical  flow  processes  reach  the  test  section,  preventing  target  flow  conditions  from  being  achieved. 
Table  2  summarises  a  Mach  13  high  total  pressure  flow  condition  which  was  attempted  with  X2.  Figure  20 
shows  several  shock  speeds  measured  at  different  points  along  the  tunnel  and  compares  these  to  theoretical 
estimates  based  on  classical  analytical  calculations  (i.e.  the  original  target  shock  speeds). 


Symbol 

Value 

Units 

Description 

PA,  0 

1.1 

MPa 

Reservoir  fill  pressure,  Air. 

Pd,  o 

30.0 

kPa 

Driver  fill  pressure,  100%  Helium. 

Prupt 

15.5 

MPa 

Primary  diaphragm  rupture  pressure. 

A 

42.5 

f-f 

Driver  compression  ratio. 

mv 

35.0 

kg 

Piston  mass. 

Psec 

150 

kPa 

Secondary  driver  fill  pressure.  Helium. 

Pshk 

330 

kPa 

Shock  tube  fill  pressure,  Air. 

Pace 

254 

Pa 

Acceleration  tube  fill  pressure,  Air. 

Lsec 

3.424 

m 

Secondary  driver  tube  length. 

Lshk 

1.301 

m 

Shock  tube  length. 

Lacc 

4.254 

m 

Acceleration  tube  length. 

M 

13.4 

- 

Predicted  Mach  number  at  nozzle  exit;  target  =  13.0. 

u 

3.950 

km/s 

Predicted  flow  velocity  at  nozzle  exit;  target  =  3.952  km/s. 

Po 

1,450 

MPa 

Predicted  total  pressure  at  nozzle  exit. 

ttt 

0.063 

ms 

Predicted  test  time. 

Table  2:  Mach  13  calculated  flow  condition. 
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Figure  20:  X2  shock  speeds  for  Mach  13  condition,  using  35  kg  piston  with  100%  Helium  driver.  Data  points 
denoted  ‘X2s...’  indicate  experimentally  measured  shock  speeds. 


Referring  to  Figure  20,  shock  speeds  are  seen  to  significantly  attenuate,  particularly  in  the  high  pressure 
(330  kPa  initial  fill  pressure)  air- filled  shock  tube.  An  analysis  from  Lld2,  which  shows  predicted  shock 
speed  as  the  incident  shock  traverses  the  length  of  the  tunnel,  demonstrates  good  agreement  with  experimen¬ 
tal  results.  This  shock  attenuation  results  in  significantly  reduced  speed  and  total  pressure  in  the  test  gas 
compared  to  the  target  flow  condition.  To  address  this  problem,  a  new  lightweight  piston  was  developed  for 
use  in  a  tuned  driver  configuration,  described  in  the  next  sections. 

It  is  noted  that  a  secondary  driver  was  used  for  this  experiment.  The  secondary  driver  is  an  intermediate 
tube  filled  with  Helium,  located  between  the  primary  diaphragm  and  the  air-filled  shock  tube.  It  is  configured 
so  that  the  sound  speed  of  the  shock  processed  Helium  in  the  secondary  driver  exceeds  the  sound  speed  of 
the  expanded  driver  gas.  This  increase  in  sound  speed  across  the  interface  between  the  two  gases  prevents 
transmission  of  transverse  acoustic  noise  in  the  driver  gas  into  the  adjacent  gas.  If  the  air  filled  shock  tube 
is  used  directly  adjacent  to  the  primary  diaphragm,  this  sound  speed  increase  is  not  achieved,  resulting  in 
significant  unsteadiness  in  the  test  flow.  This  phenomenon  is  detailed  in  Pauli  and  Stalker  [45]. 

4.3.2  Tuned  Piston  Operation 

The  concept  of  tuned  piston  operation  was  originally  proposed  by  Stalker  in  [44]  and  [46]  and  attempts  to 
increase  the  duration  over  which  driver  gas  is  maintained  at  a  useful  pressure.  It  involves  configuring  the 
driver  so  that  diaphragm  rupture  occurs  while  the  piston  still  has  sufficient  velocity  to  compensate  for  driver 
gas  loss  to  the  shock  tube  [44].  Ignoring  wave  processes  in  the  driver,  there  is  a  reference  piston  speed, 
Uref,  which  will  exactly  compensate  for  driver  gas  loss  into  the  shock  tube,  thus  resulting  in  approximately 
constant  pressure  in  the  driver.  The  actual  piston  speed  at  the  moment  of  diaphragm  rupture,  urupt,  is  non- 
dimensionalised  by  this  reference  speed,  Uref,  to  produce  Itoh’s  [47]  piston  over-drive  parameter,  \i\ 

/?  =  ^  (126) 


Stalker  [46]  proposed  the  idea  of  configuring  the  driver  such  that  /3  >  1,  thereby  “over-driving”  the 
piston.  For  /3  >  1,  the  piston  will  actually  momentarily  continue  to  increase  the  driver  pressure  following 
diaphragm  rupture,  before  pressure  begins  to  fall  again.  The  duration  of  time  over  which  this  variation  in 
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driver  pressure  is  within  acceptable  limits  (typically  considered  to  be  around  10%  of  the  target  pressure 
[48,  46,  44,  47]),  can  correspond  to  a  significantly  extended  period  of  useful  supply  time.  This  concept  is 
explained  schematically  in  Figure  21. 


Figure  21:  Effect  of  piston  over-driving  on  driver  pressure. 


4.3.3  Piston  Soft  Landing  Condition 

Over-driving  the  piston  results  in  the  piston  having  a  relatively  large  velocity  (typically  100  —  300  m/s)  at 
the  moment  when  the  diaphragm  ruptures.  Flowever,  it  is  also  necessary  to  stop  the  piston  before  it  collides 
with  the  end  of  the  compression  tube,  which  can  prove  challenging  since  the  distance  available  to  decelerate 
the  piston  is  relatively  small  for  the  high  compression  ratios  required  for  driver  performance.  Itoh  et  al. 
[47]  identified  the  types  of  motion  possible,  after  diaphragm  rupture,  as  the  piston  approaches  the  end  of 
the  compression  tube.  These  are  shown  in  Figure  22  and  are  defined  as  being  either  ‘piston  rebound’,  ‘soft 
landing’,  or  ‘direct  impact’.  The  eventual  piston  motion  depends  primarily  on  the  properties  and  initial  fill 
pressures  of  the  reservoir  and  driver  gases,  the  piston  mass,  and  the  geometry  of  the  compression  tube  and 
reservoir.  Itoh  [47]  proposes  targeting  the  soft  landing  condition  and  sizing  the  piston  buffer  so  that  it  catches 
the  piston  at  its  inflection  point  (where  piston  velocity  and  acceleration  are  simultaneously  zero;  i.e.  um  =  0 
per  Figure  22). 

A  soft  landing  condition  was  targeted  for  the  new  X2  free  piston  driver.  It  was  considered  impractical 
to  incorporate  brakes  in  the  piston  (which  help  prevent  the  rebound  motion  identified  in  Figure  22),  and 
survivable  direct  impact  is  never  feasible  for  anything  other  than  low  speed  impacts.  An  analysis  in  accor¬ 
dance  with  Stalker  [46]  indicated  that  it  was  necessary  to  make  the  new  piston  as  light  as  possible.  Structural 
strength  and  facility  interface  requirements  (i.e.  the  ability  to  use  the  piston  with  the  existing  compression 
tube  and  with  the  existing  launcher  arrangement)  placed  restrictions  on  how  light  the  piston  could  be  made. 
However,  the  final  mass  of  10.5  kg  was  determined  to  be  sufficiently  low  to  achieve  a  tuned  driver  condition 
which  would  have  sufficient  performance  to  achieve  the  target  flow  conditions.  The  new  lightweight  10.5  kg 
piston  is  shown  in  Figure  23. 

Considered  qualitatively,  tuned  free -piston  driver  conditions  require  comparatively  light  pistons  for  the 
following  reasons: 

1 .  At  the  point  of  diaphragm  rupture,  the  piston  velocity  needs  to  be  high  in  order  to  match  the  mass 
flow  of  driver  gas  into  the  shock  tube.  Especially  for  a  relatively  short  compression  tube  like  X2’s 
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(a)  Rebound  impact 


(b)  Soft  landing 

Um  0 


(c)  Direct  impact 

um  0 


Figure  22:  Characteristics  of  piston  motion  (taken  from  Itoh  [47]).  Subscript  m  refers  to  the  instant  when 
piston  acceleration  is  zero;  if  um  =  0  then  the  piston  has  an  inflection  point  where  it  can  theoretically  be 
‘caught’  by  an  appropriately  sized  buffer,  thus  avoiding  impact. 


(a)  Piston  body,  machined  from  7075-T6 
aluminium.  Note:  material  removal  across 
piston  skirt,  and  circumferential  pocketing. 


(b)  Final  piston  assembly.  Note:  nylon  (c)  Final  piston  assembly,  viewed  from 

chevron  seal  (yellow)  and  wear  bands  (dark  behind.  Note:  internal  cavity  which 
green),  and  copper  alloy  seal  support  ring.  interfaces  with  piston  launcher. 


Figure  23:  New  lightweight  piston  for  X2. 


9-46 


RTO-EN-AVT-186 


CFD  Tools  for  Design  and  Simulation 
of  Transient  Flows  in  Hypersonic  Facilities 


(approximately  4.5m),  a  piston  must  be  light  enough  to  accelerate  to  a  high  speed  over  this  short 
distance,  or  else  reservoir  pressures  must  become  prohibitively  high. 

2.  For  large  compression  ratios,  the  distance  between  the  piston  and  the  end  of  the  compression  tube  is 
quite  short.  For  a  given  driver  gas  pressure  at  rupture,  the  piston  needs  to  be  very  light  to  decelerate  to 
rest  over  this  short  distance. 

4.3.4  Calculation  of  New  Free-Piston  Driver  Conditions 

There  are  practically  limitless  combinations  of  parameters  which  will  lead  to  tuned  operation  of  a  free  piston 
driver,  but  several  design  constraints  reduce  the  design  space  to  a  more  manageable  scale: 

1.  Piston  mass:  minimum  piston  mass  is  limited  by  structural  and  interface  requirements  (10.5  kg  for 
X2’s  new  lightweight  piston). 

2.  Driver  pressure:  the  compression  tube  is  limited  by  the  magnitude  of  pressure  it  can  structurally  con¬ 
tain  (40  MPa  for  X2). 

3.  Reservoir  pressure:  the  reservoir  fill  pressure,  which  accelerates  the  piston  down  the  compression  tube, 
is  limited  by  reservoir  structural  strength  (X2  has  recently  been  temporarily  rated  to  8  MPa  to  permit 
operation  of  these  driver  conditions,  however  it  has  been  designed  for  10  MPa  and  will  be  re -rated 
accordingly  at  a  future  date). 

4.  Compression  tube  length  and  diameter:  there  is  significant  expense  involved  with  changing  the  funda¬ 
mental  configuration  of  the  facility,  therefore  compression  tube  geometry  was  assumed  to  be  fixed. 

Several  variables  remain  available  for  driver  condition  design: 

1.  Reservoir  fill  pressure  (0-8  MPa). 

2.  Driver  fill  pressure  (<1  MPa). 

3.  Driver  gas  composition  (Helium  and  Argon).  The  required  piston  speed  for  tuned  operation  depends 
on  the  speed  of  sound  of  the  compressed  driver  gas.  Reducing  the  sound  speed  (through  the  addition 
of  Argon  to  Helium),  reduces  the  required  piston  speed,  however  shock  strength  is  also  reduced. 

4.  Primary  diaphragm  thickness  and  material  (in  this  study,  diaphragm  thickness  was  limited  to  1.2,  2.0 
and  2.5  mm  thick,  cold-rolled  steel  sheet;  each  was  pre-scored  to  0.2mm  depth;  rupture  pressures  were 
assumed  to  be  15.5,  27.9  and  35.7  MPa  respectively,  based  on  previous  testing). 

5.  Buffer  length  (the  distance  from  the  extreme  end  of  the  tube  where  the  piston  makes  contact  with  the 
buffer). 

The  process  used  to  develop  new  driver  conditions  is  outlined  in  Figure  24.  The  first  step  was  to  develop  a 
rapidly  solved  0-D  perfect  gas  analytical  model  of  the  free-piston  compression  process.  The  piston  equations 
of  motion  were  obtained  from  Hornung  [49]  and  used  to  predict  piston  motion  and  driver  pressure  before  and 
after  diaphragm  rupture  (assumed  physical  processes  are  shown  in  Figure  25).  The  0-D  model  was  used  to 
manually  identify  a  range  of  potential  tuned  driver  solutions.  The  computational  time  was  sufficiently  small 
that  each  solution  could  be  quickly  identified. 

Whilst  the  0-D  model  proved  capable  of  modeling  the  driver  compression  process  fairly  effectively,  it 
could  not  make  accurate  predictions  of  required  reservoir  gas  fill  pressure.  The  reservoir  gas  expansion 
process  was  assumed  to  be  an  ideal  unsteady  expansion  as  shown  in  Figure  25b.  With  X2,  reservoir  gas  must 
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Figure  24:  X2  free-piston  driver  condition  development  process. 


pass  through  an  area  change  and  also  through  a  slotted  launcher  (refer  Figure  26).  This  convoluted  flow  path 
has  the  effect  of  throttling  the  expansion  process,  significantly  reducing  the  strength  of  the  reservoir  pressure 
force  eventually  acting  on  the  piston.  Further,  X2’s  reservoir  has  finite  length,  and  the  unsteady  expansion 
through  the  reservoir  eventually  reflects  from  the  extreme  end  and  causes  a  further  pressure  drop.  Both  of 
these  factors  necessitate  a  much  better  predictive  tool  for  the  reservoir  gas  flow. 

Lld2  was  used  to  fine  tune  the  free-piston  driver  configuration  prior  to  any  experimental  testing.  The 
code  is  capable  of  capturing  the  longitudinal  unsteady  wave  processes  which  occur  during  piston  operation 
and  includes  piston  friction,  flow  chemistry,  and  pipe-flow  viscous  effects  along  the  tube  walls.  Gradual  area 
changes  can  be  handled  by  the  code,  however  3-D  physical  processes,  such  as  flow  through  the  launcher, 
cannot  be  directly  modeled.  To  simulate  the  effect  of  these  complex  flow  paths,  Lld2  uses  loss  regions, 
which  apply  a  loss  factor  over  a  finite  length  of  the  tube  where  an  area  contraction  etc.  is  present.  Repre¬ 
sentative  loss  factors  can  only  be  determined  from  experimental  data,  therefore  development  of  loss  factors 
must  occur  in  conjunction  with  experimental  testing.  There  is  no  guarantee  that  a  loss  region  will  model  a 
disturbed  flow  region  with  useful  accuracy;  however,  anecdotal  experience  indicates  that  the  modeling  tool 
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unsteady 


expansion 

7a,  Pa.o 

IL  > 

,  ,,  ID,  PD 

~*u  aD 

“A,0 

7a,  Pa,  cia  j 

(b)  After  piston  launch,  before  diaphragm  rupture. 
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(c)  After  diaphragm  rupture. 
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Figure  25:  Assumed  free  piston  driver  states,  pre-  and  post-  diaphragm  rupture. 


Figure  26:  Piston  launcher  for  X2  (shown  detached  from  tunnel).  Note:  the  launcher  inserts  into  piston; 
reservoir  gas  must  channel  through  the  slots  in  the  launcher,  with  significant  resultant  losses  to  the  flow. 
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is  quite  effective  once  tuned  for  a  given  test  condition.  The  Lld2  driver  geometry  used  to  model  X2  with 
the  lightweight  piston  is  shown  in  Figure  27. 
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Figure  27:  Lld2  geometric  representation  of  X2  impulse  facility.  Note:  longitudinal  scale  has  been  com¬ 
pressed  to  fit  diagram  onto  page. 
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Considering  the  existing  35  kg  free  piston  driver  for  X2,  the  driver  has  been  configured  such  that  the 
amount  of  reservoir  gas  energy  imparted  to  the  piston  is  only  a  little  greater  than  that  required  to  rupture  the 
steel  diaphragm;  for  example,  if  the  reservoir  fill  pressure  is  lowered  by  10-20%,  the  piston  will  not  have 
enough  energy  to  raise  driver  pressure  to  the  diaphragm  rupture  pressure,  therefore  the  diaphragm  will  not 
rupture.  The  result  is  that  the  piston  does  not  have  significant  energy  following  diaphragm  rupture.  Further, 
since  the  piston  is  heavy,  this  energy  is  not  associated  with  a  high  velocity,  therefore  the  risk  of  significant 
impact  velocities  into  the  end  of  the  compression  tube  are  low. 

A  key  characteristic  which  differentiates  tuned  free-piston  driver  operation  with  the  lightweight  piston 
is  that  the  piston  is  given  significantly  greater  energy  than  that  which  is  required  to  break  the  diaphragm, 
since  it  must  also  have  sufficient  energy  to  continue  to  push  driver  gas  through  the  throat  of  the  driver,  at  full 
pressure,  after  the  diaphragm  has  broken.  The  lightweight  tuned  piston  has  to  be  accelerated  to  much  higher 
velocities,  be  decelerated  over  a  very  short  distance,  and  has  significantly  greater  energy  than  that  required 
to  rupture  the  diaphragm.  The  risk  of  facility  damage  due  to  uncertainties  in  the  analysis  are  much  greater, 
therefore  predictive  tools  must  be  as  accurate  as  possible.  To  achieve  this  accuracy  with  Lld2,  a  series  of 
blanked  off  tests  was  performed. 

4.3.5  Blanked-Off  Driver  Tests 

A  blanked-off  driver  test  involves  operating  a  free-piston  driver  condition  using  a  stiff,  non-rupturing  di¬ 
aphragm,  typically  manufactured  from  thick  steel.  For  this  commissioning  process,  a  PCB  pressure  trans¬ 
ducer  was  located  in  the  diaphragm,  so  that  driver  pressure  could  be  measured  during  the  piston  compression 
process.  During  a  blanked  off  test  the  piston  bounces  back  and  forth  until  the  piston  comes  to  rest.  So  long 
as  the  driver  pressure  does  not  exceed  the  facility  pressure  limit,  no  damage  will  be  done  to  the  piston.  A 
corresponding  analysis  can  be  performed  with  Lld2.  The  Lld2  model  is  then  tuned  until  an  acceptable 
level  of  correlation  is  obtained  between  the  experimental  and  numerical  pressure  traces. 

This  methodology  is  very  effective,  since  it  allows  full  correlation  of  the  driver  pressure  trace  right  up 
until  the  moment  when  the  diaphragm  rupture  pressure  is  reached.  At  this  point  with  a  normal  experiment, 
the  diaphragm  would  then  rupture,  initiating  shock  tube  flow.  If  strong  agreement  can  be  obtained  with  the 
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blanked  off  tests,  then  it  increases  confidence  that  the  post-diaphragm  rupture  piston  dynamics  will  also  be 
predicted  with  good  accuracy. 

A  broad  analysis  of  different  driver  conditions  using  the  0-D  analytical  model,  followed  by  detailed 
analysis  with  Lld2,  eventually  led  to  three  driver  configurations  which  were  considered  to  be  feasible. 
The  three  conditions  each  used  an  80%  Helium  /  20%  Argon  driver  gas  mix.  The  difference  between  the 
conditions  was  the  thickness  of  the  cold-rolled  steel  diaphragm  for  each;  1.2  mm,  2.0  mm  and  2.5  mm.  Table 
3  details  the  three  new  driver  conditions. 


Fill  pressures 

Driver  condition  ID 

Piston  mass 

Diaphragm  thickness1 

P-4,0 

PD,He,0 

PD,Ar,0 

Buffer  Length" 

[-] 

[kg] 

[mm] 

[MPa] 

[kPa] 

[kPa] 

[mm] 

LWP- 1 .2mm-Rev-0 

10.524 

1.2 

4.94 

88.2 

22.1 

100 

LWP-2.0mm-Rev-0 

10.524 

2.0 

6.85 

74.3 

18.5 

45 

LWP-2.5mm-Rev-0 

10.524 

2.5 

6.08 

61.7 

15.5 

45 

1  Diaphragms  are  manufactured  from  cold-rolled  steel  and  pre-scored  to  0.2  mm  depth. 

2  Buffer  is  comprised  of  6x50  mm  diameter  nylon  studs. 


Table  3:  X2  lightweight  piston  finalised  driver  conditions. 


Blanked-off  tests  were  performed  for  each  condition  prior  to  performing  diaphragm  rupturing  experi¬ 
ments.  Figures  28(a-c)  compare  pressure  traces  between  Lld2  predictions  and  experimental  measurement. 
Close  correlation  is  observed  for  the  average  pressure  magnitudes.  There  is  some  difference  in  the  unsteady 
behaviour  (waviness);  it  was  found  that  Lld2  had  difficulty  predicting  the  detailed  unsteady  behaviour  of 
the  driver  pressure  through  the  shaip  area  change  to  the  primary  diaphragm.  The  Lld2  pressure  traces  are 
taken  just  before  the  compression  tube  area  reduces.  It  was  found  that  loss  factors  had  to  be  increased  from 
0.5  (which  is  used  with  the  existing  35  kg  piston  Lld2  model)  to  approximately  3.5  for  the  lightweight 
piston,  to  obtain  good  agreement  between  numerical  and  experimental  driver  pressure  traces.  This  is  not 
surprising,  since  the  reservoir  pressures  are  almost  an  order  of  magnitude  higher,  and  the  piston  velocity  and 
acceleration  are  also  much  higher. 

It  is  also  noted  that  for  blanked-off  experimental  tests  with  2.0  and  2.5  mm  steel  diaphragms,  the  driver 
pressure  was  scaled  upwards  to  ensure  peak  pressure  did  not  exceed  the  facility  limit  of  40  MPa.  Since 
reservoir  pressure  has  proven  most  difficult  to  predict  accurately,  the  reservoir  pressure  was  not  scaled.  Prior 
to  the  rapid  increase  in  driver  pressure  as  the  piston  nears  the  end  of  its  stroke,  piston  dynamics  is  primarily 
dependent  on  reservoir  pressure  (i.e.  driver  pressures  are  low  for  most  of  the  piston  stroke).  Therefore  these 
scaled  blanked-off  tests  still  permit  reasonable  verification  of  most  of  the  compression  process. 

Since  Lld2  uses  a  pipe  flow  model  to  calculate  heat  loss,  it  does  not  predict  heat  loss  well  for  a  compres¬ 
sion  process  where  the  gas  is  very  hot,  but  only  moving  with  relatively  low  velocity  (i.e.  for  a  heavy,  slow 
piston).  For  these  experiments,  the  volumetric  compression  ratio  of  the  driver  gas  at  the  end  of  the  piston 
stroke  was  measured  experimentally  using  sacrificial  soft  metal  rods  fixed  into  the  end  of  the  tube,  and  it  was 
found  that  the  volumetric  compression  ratio  was  well  approximated  by  the  L 1  d2  simulations.  This  indicated 
that  heat  loss  was  not  significant  during  the  compression  process. 

Figure  29  shows  the  Lld2  predicted  piston  velocity-displacement  trajectory  for  driver  condition  LWP- 
2.0mm-Rev-0  from  Table  3.  It  can  be  seen  that  the  deceleration  of  the  piston  prior  to  reaching  the  inflection 
point  is  significant  and  that  incorrectly  locating  the  buffer  too  far  forward  of  the  tube  end  may  result  in  very 
high  speed  impact.  Driver  heat  loss  is  very  important  in  this  respect,  since  significant  heat  loss  will  result 
in  a  smaller  driver  gas  volume  at  high  pressure  and,  if  not  properly  modeled,  may  result  in  the  buffer  being 
located  too  far  forward. 
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(a)  Driver  condition  LWP-1.2mm-Rev-0  (100%  driver  fill  pressure  /  100%  reservoir  fill  pressure) 


(b)  Driver  condition  LWP-2.0mm-Rev-0  (162%  driver  fill  pressure  /  100%  reservoir  fill  pressure) 


(c)  Driver  condition  LWP-2.5mm-Rev-0  (170%  driver  fill  pressure  /  100%  reservoir  fill  pressure) 


Figure  28:  Comparison  of  experimental  and  numerical  driver  pressures  for  new  tuned  lightweight  piston 
driver  conditions  (refer  Table  3).  Experimental  pressure  traces  have  been  time-shifted  to  match  Lld2  pre¬ 
dictions. 
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x[m] 


Figure  29:  Piston  velocity  vs.  position  along  compression  tube,  driver  condition  LWP-2.0mm-Rev-0.  To 
achieve  a  soft  landing  for  the  above  example,  a  buffer  should  be  located  at  the  inflection  point  ( x  =  4.5  m). 
This  is  the  position  where  the  piston  temporarilly  comes  to  rest  before  being  pushed  forwards  again  by 
residual  reservoir  gas  pressure. 

4.3.6  Shock  Speeds  with  the  Tuned  Driver 

Figures  30(a-c)  show  example  experimental  shock  speeds  for  each  of  the  three  driver  conditions  described  in 
Table  3.  It  can  be  seen  that  there  is  no  longer  the  characteristic  shock  attenuation  observed  with  the  previous 
driver  (refer  Figure  20).  With  the  2.0mm  and  2.5mm  thick  diaphragm  conditions,  target  shock  speeds  are 
approached,  thus  achieving  the  original  goals  of  the  study. 

It  is  noted  that  the  three  driver  conditions  detailed  in  Table  3  were  achieved  without  damage  to  the 
facility.  Nylon  rods  were  used  as  the  buffer  to  catch  the  piston;  they  are  easily  cut  to  a  suitable  length 
and  have  a  high  energy  absorbing  capacity.  None  of  these  new  driver  conditions  caused  damage  to  the  nylon 
rods,  indicating  that  the  combined  analytical/numerical/experimental  development  process  managed  to  safely 
determine  tuned,  workable,  driver  conditions.  This  case  study  on  the  commissioning  of  the  new  lightweight 
piston  for  X2  has  demonstrated  a  process  which  can  be  used  to  safely  develop  new  driver  configurations,  and 
emphasises  the  use  of  the  simulation  code  L 1  d2  to  fine  tune  piston  response  prior  to  experimentation. 
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(a)  Driver  condition  LWP-1.2mm-Rev-0 


x-location  (m) 

(b)  Driver  condition  LWP-2.0mm-Rev-0 


x-location  (m) 

(c)  Driver  condition  LWP-2.5mm-Rev-0 


Figure  30:  Comparison  of  experimental,  analytical  and  numerical  shock  speeds  for  new  tuned  lightweight 
piston  driver  conditions  (refer  Table  3).  Data  points  denoted  ‘X2s...’  indicate  experimentally  measured  shock 
speeds. 
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4.4  Hybrid  simulation  of  the  X2  expansion  tube 

As  discussed  in  Section  2.2,  there  are  three  main  processes  that  occur  in  an  expansion  tube: 

1 .  free-piston  compression  and  primary  diaphragm  rupture, 

2.  high  pressure,  low  Mach  number  shock  tube  flow  and  secondary  diaphragm  rupture,  and 

3.  low  pressure,  high  Mach  number  acceleration  tube,  hypersonic  nozzle  and  dumptank  expansion. 

For  superorbital  re-entry  conditions,  the  first  two  stages  involve  all  the  important  flow  processes  occurring 
at  reasonably  high  pressures  (p  >  100  kPa)  and  low  Mach  numbers  (M  <  4),  while  the  expansion  of  the 
test  gas  into  the  acceleration  tube,  nozzle  and  dump-tank  quickly  results  in  low  pressures  (p  ~  lOkPa)  and 
high  Mach  numbers  (M  ~  10).  The  hybrid  simulation  strategy  proposed  by  Jacobs  et  al  [50]  makes  use  of 
these  flow  characteristics  to  simulate  the  X2  expansion  tube  in  an  efficient  yet  accurate  fashion.  The  in-house 
quasi-one-dimensional  Lagrangian  code  Lld2  (see  Section  3.1)  performs  well  in  the  high  pressure  and  low 
Mach  number  regimes  where  radial  flow  variation  is  small,  and  is  therefore  used  for  simulating  the  facility  up 
to  and  including  the  secondary  diaphragm  rupture.  As  the  test  gas  undergoes  the  unsteady  expansion  into  the 
acceleration  tube  post  secondary  diaphragm  rupture,  the  pressure  drops  drastically  and  the  flow  accelerates. 
The  assumption  of  a  radially  uniform  flow  subsequently  breaks  down  as  the  boundary  layer  development 
limits  the  shock  propagation  in  accordance  with  the  theory  of  Mirels  [51].  The  axisymmetric  Navier-Stokes 
equations  are  therefore  solved  with  the  mbcns2  code  (see  Section  3.2)  to  simulate  this  critical  region  of  the 
facility.  This  hybrid  simulation  technique  has  been  applied  with  success  to  an  Earth  re-entry  condition  in 
the  X3  facility  [50]  and  a  Titan  aerocapture  condition  in  the  X2  facility  [52].  Here  the  hybrid  simulation 
technique  is  applied  to  a  25  MJ/kg  CO2-N2  expansion  tunnel  condition,  with  additional  thermochemical 
nonequilibrium  analyses  of  the  secondary  diaphragm  rupture  and  steady  nozzle  expansion. 

4.4.1  Operating  Condition 

A  summary  of  the  initial  fill  conditions,  experimentally  measured  shock  speeds  and  computationally  simu¬ 
lated  freestream  conditions  are  shown  in  Table  4. 

The  Mars  test  gas  was  conservatively  taken  to  be  96%  CO2  -  4%  N2  by  volume,  as  recommended 
for  the  ESA  Radiation  Workshop  test  case  TC2-M1  [53].  A  sample  population  of  23  shots  with  Pitot  and 
static  pressure  measurements  of  the  test  flow  are  taken  to  be  representative  of  the  nominal  condition  —  the 
experimental  values  shown  in  Table  4  are  the  means  of  this  population.  The  calculated  freestream  conditions 
are  obtained  from  the  hybrid  simulation  technique  to  be  described  in  the  Section  4.4.2. 

4.4.2  Formulation  of  Hybrid  Simulation 

The  computational  approach  for  expansion  tunnel  simulation  implemented  here  considers  two  distinct  stages: 

1 .  One-dimensional  Lagrangian  simulation  of  the  high  pressure,  low  Mach  number  shock  tube  flow  and 
secondary  diaphragm  rupture,  and 

2.  Axisymmetric  Navier-Stokes  simulation  of  the  low  pressure,  high  Mach  number  acceleration  tube, 
hypersonic  nozzle  and  dumptank  expansion. 

The  driver  gas  conditions  at  primary  diaphragm  rupture  are  estimated  via  an  idealised  analysis  and  there¬ 
fore  the  piston  compression  is  not  simulated.  In  addition  to  the  main  flow  simulation,  separate  analyses 
of  thermochemical  nonequilibrium  during  the  secondary  diaphragm  rupture  and  steady  nozzle  expansion 
are  performed.  The  regions  of  the  X2  facility  considered  in  each  simulation  component  are  illustrated  in 
Figure  3 1 . 
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Experimental  parameters 

Hybrid  CFD  results 

Fill  conditions 

Reservoir 

1.28  MPa  Air 

Simulated  shock  speeds 

Primary  shock,  US3t 

3300  ±50  m/s 

Compression  tube 

30kPa  25.0%  Ar-75.0%  He 

Secondary  shock,  Usat 

6310  ±  220  m/s 

Shock  tube 

3.5  kPa  96%  C02-4%  N2 

Simulated  freestream  conditions  t 

Acceleration  tube 

9  Pa  Air 

Total  enthalpy,  htotal 

24.7  ±  0.3  MJ/kg 

Diaphragms 

Velocity,  u 

6400  ±  50  m/s 

Primary  diaphragm 

Scored  1.2  mm  steel 

Pitot  pressure,  ppitot 

65  ±  lOkPa 

Secondary  diaphragm 

13  [im  Mylar  (or)  10  pm  A1 

Static  pressure,  pstatic 

350  ±50  Pa 

Measured  shock  speeds 

Density,  p 

1.6  ±0. 10g/m3 

Primary  shock,  US)St 

3240  ±50  m/s 

Transrotational  temp.,  Ttr 

800  ±  100  K 

Secondary  shock,  Uspt 

6340  ±  230  m/s 

Vibroelectronic  temp.,  Tve 

1050  ±  100  K 

Measured  freesteam  conditions 

Equivalent  equil.  temp.,  T 

850  ±  100  K 

Pitot  pressure,  ppitot 

85  ±20kPa 

CO2  mole  fraction,  Xco2 

0.36  ±  0.04 

Static  pressure  *,  pstatic 

500  ±  170  Pa 

Mach  number,  M 

12.5  ±0.5 

Test  time,  ttest 

150  ps 

Unit  Reynolds  Number,  Re/L 

2.7  xlO4  m 1 

Table  4:  Fill  conditions,  shock  speeds  and  freestream  conditions  for  the  25  MJ/kg  CO2-N2  expansion  tunnel 
condition  in  the  X2  facility. 

Experimental  freestream  static  pressure  measured  at  wall  of  nozzle  exit. 
t  Simulated  freestream  conditions  are  averaged  over  the  central  100  mm  diameter  core  flow  and  first  100  [is 
of  test  time. 

4.5  Free-piston  compression  and  shock  tube  flow 

The  quasi-one-dimensional  nature  of  the  L 1  d2  geometry  can  result  in  boundary  layer  heat  loss  being  sig¬ 
nificantly  under-estimated.  Previous  attempts  at  modeling  the  free-piston  compression  of  the  driver  gas  have 
demonstrated  this  inadequacy  [54],  with  the  driver  gas  temperature  at  primary  diaphragm  rupture  being  un¬ 
reasonably  high.  The  Lld2  simulation  of  the  shock  tube  flow  is  therefore  begun  at  the  moment  of  primary 
diaphragm  rupture.  The  driver  gas  pressure  and  slug-length  at  rupture  are  obtained  from  an  idealised  model 
of  the  free-piston  compression,  whilst  the  temperature  is  obtained  parametrically  by  matching  the  experi¬ 
mentally  measured  primary  shock  speed.  Momentum  loss  at  the  primary  diaphragm  station  area  change  is 
accounted  for  through  a  loss-per-unit-length  factor  Kl  of  0.25.  The  CO2-N2  test  gas  is  described  by  an 
equilibrium  equation-of-state  using  the  curve  fits  provided  by  the  CEA  program  [55],  as  are  the  respective 
transport  coefficients.  The  Ar-He  driver  gas  is  described  as  a  mixture  of  ideal  gases. 

4.6  Light  secondary  diaphragm  rupture 

When  the  primary  shock  through  the  stagnant  test  gas  reaches  the  light  secondary  diaphragm,  the  first  few 
millimeters  of  shock-processed  test  gas  are  stagnated  by  the  resulting  reflected  shock.  For  the  condition  at 
hand  the  ratio  of  the  stagnated  to  freestream  density  of  the  test  gas  is  approximately  104  :  1  —  taking  area 
change  into  account  and  using  the  freestream  conditions  from  Table  4,  the  observed  150  [is  of  test  flow 
originates  from  within  the  first  2  mm  of  stagnated  test  gas.  It  follows  that  the  subsequent  expansion  of  this 
very  small  test  gas  volume  through  the  acceleration  tube  and  nozzle  must  be  modeled  accurately. 

Diaphragms  are  typically  modeled  in  Lld2  as  fixed-wall  boundary  conditions  that  are  held  in  place 
for  a  small  period  of  time  to  allow  the  reflected  shock  to  form  before  being  released.  This  method  is  re¬ 
ferred  to  as  the  holding-time  model,  and  is  widely  implemented  due  to  its  relative  simplicity.  Bakos  and 
Morgan  [56]  demonstrated  the  limitations  of  the  holding-time  model  through  comparison  with  an  inertial 
diaphragm  model.  Space-time  representations  of  both  the  holding-time  and  inertial  diaphragm  models  are 
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Figure  31:  Regions  of  the  X2  facility  modelled  in  the  present  hybrid  simulation  technique 


(a)  Holding  time  model 


Figure  32:  Space-time  diagrams  of  light  secondary  diaphragm  rupture  in  an  expansion  tube  (adapted  from 
Bakos  and  Morgan  [56]). 


shown  in  Figure  32.  The  inertial  diaphragm  model  represents  the  diaphragm  as  an  ideal  piston  —  the  di¬ 
aphragm  shears  cleanly  at  the  tube  wall,  remains  planar  during  its  motion  and  provides  only  inertial  resistance 
to  the  expanding  test  gas.  Where  a  test  gas  particle  immediately  upstream  of  the  secondary  diaphragm  ex¬ 
periences  an  infinite  rate  of  expansion  under  the  holding-time  model,  the  inertial  diaphragm  model  shows 
the  contact  surface  must  accelerate  to  an  asymptotic  limit  over  a  finite  period  of  time.  The  duration  a  test 
gas  particle  spends  in  the  unsteady  expansion  directly  determines  the  degree  of  thermochemical  relaxation 
experienced  —  high  expansion  rates  will  result  in  frozen  thermochemistry,  while  low  expansion  rates  will 
tend  towards  equilibrium. 

The  light  secondary  diaphragms  used  for  the  sample  population  shots  were  13  pm  Mylar  sheets  weighing 
approximately  70  pp  —  the  assumption  that  this  sheet  will  stay  intact  through  the  entire  test  gas  expansion 
is  clearly  invalid.  A  decaying  inertial  diaphragm  model  is  therefore  proposed  for  the  present  analysis,  where 
the  diaphragm  reduces  in  mass  exponentially  after  coming  in  contact  with  the  hot  stagnated  test  gas.  The 
expression  for  the  diaphragm  mass  rn  rate-of-change  is  taken  to  be: 

dm 
dt 


/decay  x  m  for  m  >  miimit 
0  for  m  <  miimit 


(127) 
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where  /decay  is  a  time  constant  and  rniimit  an  imposed  lower  mass  limit  to  prevent  numerical  instabilities.  Sim¬ 
ulations  of  the  secondary  diaphragm  rupture  and  subsequent  unsteady  expansion  are  performed  for  holding¬ 
time,  inertial  diaphragm  and  decaying  inertial  diaphragm  models  with  finite-rate  chemistry  using  the  Lid 
code. 


4.7  Acceleration  tube,  nozzle  and  dumptank  flow 

The  axisymmetric  flow  solver  mbcns  is  used  for  the  low  pressure,  high  Mach  number  acceleration  tube, 
nozzle  and  dumptank  flow.  The  computational  domain  is  subdivided  into  a  number  of  interconnected  blocks, 
as  illustrated  in  Figure  33,  which  allows  the  calculation  to  be  parallelised  over  a  distributed  memory  cluster 
computer. 

0.3  Boundary  conditions: 

— fixed  temperature  wall  (T  =  296  K) 

.  extrapolated  outflow 

0.2.  -  symmetry  axis 

- block  boundary 


Inflow  Acceleration  tube  Hypersonic  nozzle  Dump-tank 

plane 


[60]  ]  [61] 


Figure  33:  Computational  domain  for  axisymmetric  simulation  of  the  X2  expansion  tunnel  downstream  of 
the  secondary  diaphragm. 


Although  the  holding-time  secondary  diaphragm  rupture  model  is  expected  to  underestimate  test  gas 
recombination  through  the  unsteady  expansion,  it  provides  a  simple  pre -rupture  flowfield  that  can  be  eas¬ 
ily  used  to  initialise  the  mbcns  simulation.  Using  a  holding-time  of  10 /rs  the  stagnated  test  gas  slug  is 
calculated  to  be  4  mm  long  at  rupture.  A  4  mm  block  thus  precedes  the  acceleration  tube,  as  shown  in  Fig¬ 
ure  33,  which  is  filled  with  the  equilibrium  stagnated  test  gas  conditions.  The  upstream  face  of  this  block  is 
a  transient  inflow  boundary  condition,  applying  the  one-dimensional  flow  solution  from  the  Lid  simulation 
uniformally  over  the  inflow  plane. 

Although  thermal  nonequilibrium  is  anticipated  to  occur  during  the  nozzle  expansion  process,  the  inclu¬ 
sion  of  both  thermal  and  chemical  nonequilibrium  for  such  a  large-scale  simulation  is  computationally  pro¬ 
hibitive.  Thermal  equilibrium  is  therefore  assumed  for  the  Navier-Stokes  simulations  of  the  post-secondary 
diaphragm  rupture  flow,  while  separate  inviscid  simulations  of  the  hypersonic  nozzle  are  performed  with 
the  two-temperature  model  proposed  by  Park  [57,  58].  The  stagnated  test  gas  upstream  of  the  secondary 
diaphragm  just  prior  to  rupture  contains  negligible  ionic  species,  and  therefore  only  the  neutral  dissociation 
and  exchange  reactions  listed  in  Table  5  (reactions  1  to  16)  are  included  in  the  Navier-Stokes  simulations. 

4.7.1  Results  and  comparison  with  experiment 

Lagrangian  shock  tube  simulations  A  driver  slug  length  of  170  mm  and  piston  velocity  of  50  m/s  at  pri¬ 
mary  diaphragm  rupture  was  determined  through  an  idealised  model  of  the  free-piston  compression.  A  driver 
slug  temperature  of  2800  K  was  selected  through  the  matching  of  experimental  conditions  in  a  parametric 
analysis.  A  comparison  of  the  Lld2  and  experimental  static  pressure  history  from  shot  x2s248  is  shown  in 
Figure  34a.  Shock  arrival  time  observed  in  experiment  at  each  transducer  location  is  matched  by  the  simula¬ 
tions  within  experimental  and  modeling  uncertainties.  Although  the  experimental  data  shows  a  higher  initial 
pressure  rise  for  transducer  ST1,  excellent  agreement  is  shown  for  the  first  800  ns  of  flow  0.5m  downstream 
at  ST3.  A  space-time  representation  of  logarithmic  pressure  contours  is  shown  in  Figure  34b.  The  tail  of 
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Neutral  reactions 

Activation 

energy 

Ionic  reactions 

Activation 

energy 

Dissociation  reactions  Associative  ionisation  reactions 


1  C02  +  M  <=^  CO  +  O  +  M 

2  CO  +  M  4=>-  C  +  O  +  M 

3  N2  +  M  4=>  N  +  N  +  M 

4  02  +  M  4=>  O  +  O  +  M 

5  NO  +  M  N  +  O  +  M 

6  CN  +  M  C  +  N  +  M 

7  C2  +  M  <==>  C  +  C  +  M 
Exchange  reactions 

8  NO  +  O  4=>  N  +  02 

9  N2  +  O  4=>-  N  +  02 

10  co  +  o^c  +  o2 

11  co  +  c^c2  +  o 

12  CO  +  N^-CN  +  O 

13  N2  +  C<=^CN  +  N 

14  CN  +  O  NO  +  C 

15  CN  +  C  4=>  C2  +  N 

16  C02  +  0^02  +  C0 


-526  kJ/mol 
-1073  kJ/mol 
-941  kJ/mol 
-497  kJ/mol 
-628  kJ/mol 
-590  kJ/mol 
-581  kJ/mol 

-162  kJ/mol 
-319  kJ/mol 
-575  kJ/mol 
-48  kJ/mol 
-321  kJ/mol 
-193  kJ/mol 
-121  kJ/mol 
- 1 1  kJ/mol 
-231  kJ/mol 


17 

N  +  O  4= 

=4  NO+  +  e“ 

-265  kJ/mol 

18 

0  +  04= 

=^0+  +  e- 

-670  kJ/mol 

19 

C  +  O  4= 

=4  CO+  +  e- 

-275  kJ/mol 

Charge  exchange  reactions 

20  NO+  +  C  <=^  NO  +  C+  - 1 93  kJ/mol 

21  O^  +  O  <=^  0+  +  02  -150  kJ/mol 

22  NO+  +  N  4==>-  0+  +  N  -106  kJ/mol 

23  NO+  +  O  4=A  O^  +  N  -404  kJ/mol 

24  CO  +  C+  4=>  CO+  +  C  -261  kJ/mol 

25  02  +  C+  <=^  O^  +  C  -78  kJ/mol 

Electron-impact  ionisation  reactions 

26  C  +  e-  4=>  C+  +  e“  +  e_  -1087  kJ/mol 

27  O  +  e~  4=>-  0+  +  e_  +  e-  -1318  kJ/mol 


Table  5:  Neutral  and  ionic  chemical  reactions  from  the  C02  -  N2  reaction  scheme  of  Park  et  al  [58]  (Note 
that  the  photo-ionisation  reactions  have  been  omited). 


the  expansion  fan  eminating  from  the  primary  diaphragm  station  can  be  seen  to  reflect  off  the  compression- 
shock  tube  area  change  and  eventually  coalesce  with  the  shock  at  t  =  700  //s.  Stagnation  conditions  behind 
the  reflected  shock  at  diaphragm  rupture  were  determined  to  be  10.3  MPa  and  4370  K  with  an  equilibrium 
CO2  fraction  of  30%  by  mass. 

Secondary  diaphragm  rupture  analysis  Holding-time,  inertial  diaphragm  and  decaying  inertial  diaphragm 
models  were  compared  in  the  analysis  of  secondary  diaphragm  rupture.  Space-time  plots  of  shock  propa¬ 
gation  obtained  from  finite-rate  Lid  simulations  with  the  three  secondary  diaphragm  models  considered  are 
shown  in  Figure  35a  with  experimental  shock  arrival  times  overlaid.  The  experimental  transducer  shock  ar¬ 
rival  times  are  bounded  by  the  holding-time  and  inertial  diaphragm  solutions,  with  the  holding-time  model 
predicting  the  strongest  secondary  shock  as  expected.  Implementing  a  decaying  inertial  diaphragm  model 
with  a  time  constant  f decay  of  5  x  104  s_1  was  found  to  match  the  secondary  shock  propagation  observed  in 
experiment. 

CO2  mass  fraction  distributions  at  t  =  trupture+ 150  ps  for  each  of  the  Lagrangian  simulations  conducted 
are  shown  in  Figure  35b.  The  test  gas  composition  was  found  to  be  chemically  frozen  from  this  time  onwards 
for  all  diaphragm  models.  The  holding-time  model  gives  the  lowest  levels  of  C02  recombination,  with  a 
much  more  pronounced  variation  over  the  slug  when  compared  to  the  relatively  uniform  inertial-diaphragm 
results.  An  unrealistically  large  initial  rate  of  expansion  for  the  test  gas  particles  immediately  upstream 
of  the  secondary  diaphragm  gives  insufficient  time  for  any  recombination  to  occur  before  freezing  occurs. 
The  holding-time  recombination  levels  only  begin  to  approach  that  of  the  inertial  diaphragm  simulations  at 
the  very  end  of  the  test  slug.  The  inertial  diaphragm  model  with  no  mass  decay  shows  significantly  higher 
levels  of  CO2  recombination  and  a  more  uniform  distribution.  The  decaying  inertial  diaphragm  model  gives 
recombination  levels  bounded  by  the  previous  two  models  as  expected,  with  an  average  C02  mass  fraction 
in  the  test  slug  of  0.46.  Assuming  that  the  decaying  inertial  diaphragm  model  closely  represents  the  physical 
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(a)  Comparision  of  simulated  and  experimental  static  pressure 
traces  in  the  shock-tube 


x2s247-FR-lld05_  log-p,Pa  LlD:space— time— plot 

xl  x2  dx  4.50e+00  1.25e+01  2.00e+00  yl  y2  dy  0.00e+00  2.50e-03  5.00e-04 
vl  v2  dv  1.00e+00  7.00e+00  1.00e-01 


O 


x,m 

(b)  Space-time  diagram  of  logarithmic  pressure  contours. 


Figure  34:  Results  from  quasi-one-dimensional  Lagrangian  shock  tube  simulation  with  equilibrium  chem 
istry. 


rupture  behaviour,  implementing  a  holding-time  rupture  model  to  initialise  the  Navier-Stokes  simulations 
will  result  in  CO2  levels  being  underpredicted  by  between  25%  and  4%  over  the  test  time. 


Navier-Stokes  expansion  tunnel  simulation  Viscous  Navier-Stokes  simulations  of  the  expansion  tunnel 
flow  were  conducted  with  finite-rate  chemistry  and  a  pre -rupture  flowfield  described  by  the  Lld2  solution 
with  a  holding-time  secondary  diaphragm.  Initial  simulations  with  37  x  4500  cells  in  the  acceleration  tube 
and  a  stagnated  test  gas  slug  length  of  2mm  resulted  in  unphysical  pockets  of  hot,  low  density  gas  forming 
along  the  axisymmetric  axis.  Increasing  the  cell  distribution  to  43  x  5155  (1  mm  squares,  average)  reduced 
this  phenomena  slightly,  however  disturbances  were  still  observed  during  the  test  time.  The  test  gas  slug 
length  was  increased  to  4  mm,  corresponding  to  a  hold-time  of  10  / is ,  and  the  disturbances  no  longer  formed 
during  the  test  time.  The  stability  of  the  flow  during  the  test  time  for  this  final  simulation  is  illustrated  in 
the  contours  of  pressure  through  the  nozzle.  Figure  36.  The  freestream  conditions  from  this  simulation  are 
shown  in  Table  4. 

A  comparison  of  the  mbcns2  solution  with  experimentally  obtained  pressure  traces  is  shown  in  Fig¬ 
ure  37.  The  simulated  secondary  shock  propagation  shown  in  Figure  37a  closely  matches  that  of  shot  x2s247 
which  is  close  to  the  mean  of  the  sample  population.  Good  agreement  with  experiment  is  shown  for  both 
static  and  Pitot  pressure  at  the  nozzle  exit  for  the  first  200  ns  of  flow,  Figures  37b  and  37c.  The  simulated 
Pitot  pressure  profile  100  mm  downstream  of  the  nozzle  exit  shows  reasonable  agreement  with  that  obtained 
through  the  experimental  Pitot  pressure  survey,  Figure  37d.  For  the  central  100  mm  of  flow  a  test  article 
is  likely  to  occupy  slight  property  variation  exists,  most  noticeably  at  the  extremities  which  has  Pitot  pres¬ 
sure  10%  higher  than  at  the  central  axis.  Overall,  the  mbcns2  simulation  shows  sufficient  correlation  with 
experiment  to  proceed  with  further  validation  through  bluntbody  shock  layer  simulations. 


Nonequilibrium  nozzle  expansion  analysis  The  averaged  nozzle  entrance  conditions  from  the  viscous, 
thermal  equilibrium  mbcns  expansion  tunnel  simulation  were  used  as  initial  conditions  for  inviscid,  thermal 
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Distance  from  secondary  diaphragm,  x  (mm) 


(a)  XT  diagram  from  ST1  to  AT2. 
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(b)  Test  gas  CO2  mass  fractions  150  [is  after  rupture. 


Tes 

-gas  slug 

Inertial  diaphragms 

Holding-time  diaphragm 

Holding-time  diaphragm 
Inertial  diaphragm,  m  =  68  pg,  fdeqa(  =  0 1/s  — 
Inertial  diaphragm,  m  =  68  pg,  fdecav  =  50000  1/s  — 

Figure  35:  Comparison  of  shock  propagation  and  test  gas  recombination  for  various  light  secondary  di¬ 
aphragm  models. 
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(a)  Static  pressure  history  at  AT3,  AT4  and  AT5. 


(b)  Static  pressure  history  at  nozzle  exit. 


^rupture  A  1175/XS. 


Figure  37:  Comparison  of  computational  and  experimental  pressures  through  the  acceleration  tube  and  noz¬ 
zle. 


nonequilibrium  simulations  of  the  hypersonic  nozzle  expansion.  Test  gas  temperature  contours  through  the 
nozzle  after  800  ps  of  flow  are  shown  in  Figure  38a,  while  the  respective  centreline  temperature  profile  is 
shown  in  Figure  38b.  Both  temperatures  along  the  centreline  are  seen  to  remain  constant  until  the  Mach 
cone  arrives  250  mm  downstream  of  the  nozzle  entrance.  The  transrotational  temperature  then  begins  to 
drop  rapidly  as  the  test  gas  expands,  while  the  vibroelectronic  temperature  relaxes  at  a  much  lower  rate. 
After  800  mm  the  transrotational  temperature  along  the  centreline  is  essentially  frozen  at  800  K  while  the 
vibroelectronic  temperature  continues  to  slowly  relax.  The  flow  close  to  the  nozzle  wall  remains  in  a  state 
of  thermal  equilibrium  at  1100K  through  the  entire  length  due  to  a  much  weaker  expansion  than  at  the 
centreline.  Slightly  increased  relaxation  begins  to  occur  as  the  test  gas  expands  into  the  dump  tank,  resulting 
in  a  final  temperatures  100  mm  downstream  of  the  nozzle  exit  of  Tve  ~  1050  K  and  Ttr  ~  800  K. 

4.7.2  Quality  of  the  Solution 

A  hybrid  simulation  technique  for  the  X2  expansion  tube  incorporating  a  one-dimensional  Lagrangian  sim¬ 
ulation  of  the  shock  tube  and  secondary-diaphragm  and  an  axisymmetric  Navier-Stokes  simulation  of  the 
acceleration  tube,  nozzle  and  dump-tank  has  been  applied  to  a  25  MJ/kg  CO2-N2  condition.  Good  agree- 
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(a)  Contours  of  temperature. 


Axial  distance  from  nozzle  entrance,  x  (m) 


(b)  Centreline  temperature  profiles. 

Figure  38:  Test  gas  temperature  distributions  through  the  hypersonic  nozzle  after  800  //s  of  flow. 


ment  with  experimentally  measured  static  pressure  traces  and  shock  speeds  was  obtained  in  both  the  shock 
and  accerlation  tubes,  and  Pitot  pressure  was  also  found  to  be  accurately  predicted  in  the  test-section.  Al¬ 
though  the  application  of  a  simplified  holding-time  model  allows  the  shock  propagation  and  pressure  levels 
to  be  accurately  reproduced  in  the  Navier-Stokes  simulation,  test  gas  recombination  following  secondary 
diaphragm  rupture  was  found  to  be  under  predicted  by  as  much  as  25  %  when  compared  to  a  simulation 
considering  an  inertial  secondary  diaphragm.  Furthermore  considerable  thermal  nonequilibrium  is  shown 
to  develop  through  the  steady  nozzle  expansion,  with  Tve  and  Ttr  deviating  by  24  %  and  -6  %  respectively 
from  the  equilibrium  temperature  Teq.  Therefore  while  the  density,  pressure  and  velocity  obtained  via  the 
hybrid  simulation  technique  is  considered  to  be  accurate,  the  species  mass-fractions  and  thermal  equilibrium 
temperature  are  considered  approximate  only. 


RTO-EN-AVT-186 


9-63 


CFD  Tools  for  Design  and  Simulation 
of  Transient  Flows  in  Hypersonic  Facilities 


ORGANIZATION 


4.8  Axisymmetric  simulation  of  the  X2  shock  tube 
4.8.1  Experiment  discription 

A  comprehensive  set  of  shock  tube  experiments  with  a  representative  Titan  test  gas  (98%  N2  &  2%  CH4  by 
volume)  over  a  pressure  and  velocity  range  of  2  to  1000  Pa  and  4  to  10.5  km/s  respectively  were  conducted 
by  Brandis  [59]  in  the  X2  facility.  Following  the  work  of  Gollan  [60,  61],  a  1  Torr  5.7  km/s  condition 
was  the  focus  of  a  recent  computational  study  [62]  where  the  CN  Violet  radiation  intensity  predicted  by  a 
one-dimensional  and  an  axisymmetric  simulation  were  compared  with  experiment.  The  aim  of  the  study 
was  to  assess  the  assumption  of  one-dimensional  variation  of  properties  during  the  test-time,  as  the  optical 
line-of-sight  passes  through  the  expanding  boundary  layer  and  Mach  cone  in  the  test  section. 

Table  6  summarises  the  fill  conditions,  shock  speeds  and  spectral  measurements  performed  in  the  X2 
experiments  that  were  simulated. 


x2s522 


x2s645 


Shot  number 


86.0%  He,  14.0%  Ar 
15.1  MPa 


Driver  gas  composition  * 
Primary  diaphragm  burst  pressure 
Test  gas  composition  * 
Shock  speed,  Us  f 
Pressure,  p ^ 

Exposure  time,  tex. 
Spectrometer  range 


98%  N2,  2%  CH4 
5697  ±  54  m/s  5658  ±  53  m/s 
133  ±  0.5  Pa 
100  ns 
308-450  nm 


Table  6:  Two  X2  shots  targeting  a  1  Torr  5.7  km/s  Titan  entry  condition 

*  Gas  percentage  compositions  are  by  volume 
t  Given  shock  speed  is  that  between  transducers  AT5  and  AT7  (see 
Figure  4) 

Two  shots  x2s522  and  x2s645  were  both  performed  with  an  initial  test  gas  pressure  of  approximately 
1  Ton-  targeting  a  shock  speed  of  5.7  km/s;  the  measured  shock  speeds  matched  this  within  the  bounds  of 
measurement  uncertainty.  Although  the  pressure  of  this  condition  is  an  order  of  magnitude  higher  than  that 
of  Huygens-type  direct  entry  peak  heating  [63],  it  is  a  useful  condition  for  investigating  the  chemical  kinetics 
of  Titan  gas  closer  to  thermochemical  equilibrium. 

4.8.2  Formulation  of  Axisymmetric  Simulation 

The  axisymmetric  Navier-Stokes  simulation  technique  implemented  here  is  based  on  that  described  in 
Ref.  [60],  however  the  driver  gas  conditions  are  assumed  to  be  uniform  (T=3450K,  p=  15 . 1  MPa,  u=0m/s) 
at  diaphragm  rupture  for  the  present  work.  A  single-temperature  gas  model  is  implemented  alongside  the 
Go  keen  [64]  chemical  reaction  scheme  where  the  ionic  reactions  have  been  omitted.  This  relatively  low  or¬ 
der  thermochemistry  model  was  used  due  to  computational  resource  constraints.  The  computational  domain 
for  the  axisymmetric  Navier-Stokes  simulation  of  shot  x2s522  is  presented  in  Figure  39.  The  computational 
domain  extends  from  the  front  piston  face  to  1  m  into  the  dump-tank;  the  remaining  length  of  the  dump-tank 
is  not  included. 

Again,  the  simulations  were  performed  with  the  mbcns2  code  described  in  Section  3.2.  The  compu¬ 
tational  domain  is  decomposed  into  160  roughly  equally  sized  blocks  -  5  in  the  compression  tube  cavity, 
123  along  the  shock-tube  and  32  in  the  dump-tank.  Two  grid  resolutions  were  considered  -  (a)  20  x  2000 
cells  in  the  shock-tube  and  (b)  40  x  4000  cells  in  the  shock-tube.  The  cell  discretization  in  the  compression 
cavity  and  dump-tank  are  matched  to  the  shock  tube  on  a  cell  per  unit  area  basis.  The  higher  resolution 
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Figure  39:  X2  computational  domain  for  axisymmetric  Navier-Stokes  simulation  of  shot  x2s522 


Axial  location,  x  (m) 


Figure  40:  Comparison  of  shock-speed  from  experiments  and  simulations  for  the  1  Torr  X2  condition 


simulation  has  a  total  of  580,180  cells.  The  simulations  were  run  on  40  quad-core  Intel  Woodcrest  nodes  and 
the  40  x  4000  simulation  required  approximately  140  hours  for  completion. 

4.8.3  Results  and  comparison  with  experiment 

Figure  40  compares  the  shock-speeds  measured  from  experiment  and  those  obtained  from  the  computational 
simulations.  The  shock  speeds  from  experiment  are  determined  by  locating  the  arrival  time  of  the  shock 
at  each  of  the  8  pressure  transducers  along  the  shock  tube  and  calculating  the  time-of-flight  in  between 
each  pair.  The  spatial  uncertainty  is  1  mm  and  the  temporal  uncertainty  is  nominally  1  /is.  Although  the 
40  x  4000  cell  simulation  agrees  with  the  experimentally  measured  shock  speed  of  approximately  5.7  km/s 
at  the  tube  exit,  the  shock  propagation  through  the  shock  tube  is  quite  different.  The  experiments  show  a 
quick  initial  rise  and  then  plateau  in  shock  speed  in  the  first  2  meters  after  the  primary  diaphragm,  whilst 
the  sudden  rise  in  the  simulations  is  only  2  meters  from  the  tube  exit.  The  sudden  rise  in  shock  speed  is 
thought  to  be  due  to  the  reflected  expansion  wave  catching-up  with  and  accelerating  the  shock.  The  location 
of  the  shock/expansion  wave  interaction  is  determined  largely  by  the  initial  temperature  of  the  driver  gas 
-  for  the  present  simulations  this  was  set  to  3,450  K  which  is  approximately  the  isentropic  compression 
temperature.  If  the  driver  gas  compression  process  is  not  adiabatic  the  temperature  of  diaphragm  rupture 
may  be  substantially  less,  and  the  the  shock/expansion  wave  interaction  location  would  be  different. 
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The  simulated  and  x2s522  static  pressure  traces  at  transducers  ST1,  AT4,  AT5  and  AT7  are  compared 
in  Figures  41a  to  4 Id  respectively.  Generally  the  quantitative  agreement  with  experiment  is  quite  good. 
Of  particular  note,  however,  is  the  discrepancy  in  expansion  wave  arrival  time  at  the  ST1  transducer  -  in  the 
experiment  the  arrival  time  is  approximately  550  /rs,  while  in  simulation  it  is  600  //s.  Such  a  large  discrepancy 
so  close  to  the  primary  diaphragm  (3.5  m)  indicates  that  the  expansion  of  the  driver  gas  is  not  being  modeled 
correctly.  By  the  AT4  transducer  the  simulated  shock  is  in  agreement  with  experiment,  and  this  is  reflected 
in  the  close  agreement  in  the  pressure  histories  at  AT4,  AT5  and  AT7.  Most  critically,  the  post-shock  static 
pressure  level  of  approximately  50kPa  at  the  final  transducer  before  the  tube  exit  is  matched  almost  exactly. 

As  described  in  Ref.  [59],  the  spectral  measurements  in  the  X2  facility  considered  in  the  present  work 
are  made  as  the  shock  is  expanding  from  the  shock  tube  exit  into  the  dump-tank  (see  Figure  3).  Figure  42 
displays  CN  mass  fraction  contours  from  the  mbcns2  simulation  of  the  1  Ton-  Titan  condition  as  the  shock 
emerges  into  the  test  section  for  spectral  measurement.  The  spectral  measurements  arc  taken  at  approxi¬ 
mately  t=1720  / / s .  Comparing  Figures  42a  and  42e,  the  expansion  of  the  shock  layer  into  the  test  section  is 
seen  to  alter  the  quasi-one-dimensional  nature  of  the  flow  present  inside  the  tube. 

To  evaluate  the  assumption  that  the  emerging  shock-layer  remains  approximately  one-dimensional  until 
the  spectral  measurements  are  taken,  one-dimensional  post-shock  relaxation  simulations  with  the  poshax 
code  are  compared  with  the  axisymmetric  mbcns2  simulations.  Both  simulations  use  the  aforementioned 
1  temperature  thermal  model  with  the  reduced  Gok£en  [64]  chemistry  scheme.  Figures  43a  through  43c 
present  comparisons  of  temperature,  CN  mole  fraction  and  radiative  intensity  profiles,  respectively,  for  the  1 
Ton-  Titan  condition  in  X2.  The  experimentally  measured  intensity  profile  from  shot  x2s645  is  overlayed  on 
Figure  43c  and  has  been  scaled  to  match  the  order  of  magnitude  of  the  one  temperature  simulations.  While 
the  one-dimensional  simulation  uses  the  Rankine-Flugoniot  relations  to  calculate  the  gas  condition  directly 
behind  the  shock  (shock-slip  is  not  considered),  the  finite-volume  Navier-Stokes  simulation  captures  the 
shock  over  a  number  of  cells  resulting  in  a  diffused  shock  front.  Furthermore  the  curvature  of  the  shock  as 
it  emerges  into  the  test  section  is  only  modelled  in  the  axisymmetric  simulations.  These  effects  are  evident 
in  Figures  43a  and  43b,  where  the  change  in  temperature  and  CN  mass  fraction  across  the  shock  are  more 
gradual  in  the  Navier-Stokes  simulations  and  peak  at  lower  values.  Consequently  the  calculated  radiation 
intensity  profile  from  the  Navier-Stokes  simulation  exhibits  a  qualitative  improvement  over  that  from  the 
one-dimension  post-shock  relaxation  simulation,  Figure  43c.  The  ratio  of  peak  to  equilibrium  intensity  is 
much  closer  in  the  Navier-Stokes  simulation  and  the  rise  to  the  peak  intensity  immediately  behind  the  shock 
is  more  realistic. 

4.8.4  Quality  of  the  Solution 

The  shock  speed  and  static  pressure  levels  calculated  at  AT5  and  AT7  matched  that  measured  in  experiment, 
however  the  shock  speed  was  under-predicted  in  the  early  stages  of  shock  propagation.  A  better  repre¬ 
sentation  of  shock  propagation  would  require  an  axisymmetric  simulation  of  the  free-piston  compression 
rather  than  the  assumption  of  uniform  driver  properties  at  primary  diaphragm  rupture.  Post-processing  of 
the  Navier-Stokes  simulation  data  resulted  in  a  radiation  intensity  profile  that  gave  improved  qualitative 
agreement  with  experiment  compared  to  that  from  a  one-dimensional  post-shock  relaxation  simulation. 
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(a)  ST1,  x=7.381  m 


Time,  t  ((is) 


(b)  AT4,  x=12.675  m 


Time,  t  (ns) 

(c)  AT5,  x=12.855  m 


Figure  41:  Comparison  of  simulated  and  experimentally  measured  static  pressure  histories  for  the  1  Torr 
condition 
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(c)  t=17 10  fis 


(d)  t=17 15  /is 


(e)  t=1720  fis 


(f)  t=1725  (is 


(g)t=1730/xs  (h)  t=1735  fj,s 


Figure  42:  CN  mass  fraction  contours  from  Navier-Stokes  simulation  of  the  1  Torr  Titan  condition  (x2s522) 
in  the  X2  facility.  The  spectral  measurements  are  taken  at  approximately  t=  1 720  //s. 
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Figure  43:  Comparison  of  one-dimensional  and  axisymmetric  (centerline,  t  =  1720  //s)  profiles  for  the  1 
Ton-  Titan  condition  in  X2 
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5  WHERE  TO  FROM  HERE? 

In  the  flow  simulation  codes  that  we  have  described  here,  the  modelling  of  unsteady  gas-dynamic  wave  mo¬ 
tion,  viscous  effects  and  finite-rate  chemical  processes  are  all  well  addressed.  Current  work  is  continuing  on 
turbulence  modelling,  radiation  energy  exchange  and  the  inclusion  of  moving  pistons  in  the  axisymmetric 
flow  simulation  code.  With  these  developments,  it  should  be  possible  (given  sufficient,  very  large  compu¬ 
tational  resources)  to  model  a  free-piston  driven  shock  tunnel  or  expansion  tube  from  end  to  end  with  an 
axisymmetric  flow  simulation. 

On  the  software  development  side,  there  is  a  new  version  of  both  Lid  and  mbcns  just  about  ready  for 
general  use.  The  thermochemistry  module  has  been  rewritten  and  generalized  so  that  it  can  be  include  a 
range  of  thermal  nonequilibrium  models  that  are  needed  to  better  handle  high-temperature  radiating  flow 
fields.  The  axisymmetric  flow  solver  has  also  been  combined  with  a  full  three-dimensional  flow  solver  to 
become  the  new  code  Eilmer3  [65]. 
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A  Nomenclature,  Units 


A 

a 

Cv 

Cp 

D 

E 

e 

F 

FP 

Fwall 

Floss 

f 

f 

H 

h 

h 

h'i 

hi 

K 

k 

L 

M 

MW 

m 

mp 

n 

n 

h,p 

P,p 

PD, He,  0 
PD,Ar,  0 

Pr 

Prupt 

Q 

Q 

q 

R 

Ro 

Re 

r 

r 

S 

St  __ 
Ul,Ur 


duct  area  or  cell  area,  m2 
speed  of  sound,  m/s 

specific  heat  at  constant  volume,  J /kg  •  K 

specific  heat  at  constant  pressure,  J /kg  •  K 

(effective)  duct  diameter,  m 

total  energy  per  unit  mass  e  +  \  u2,  J/kg 

specific  internal  energy,  J/kg 

array  of  flux  terms 

piston  friction  force,  N 

wall  shear  force  due  to  viscous  effects,  N 

effective  force  due  to  pipe  fitting  losses,  N 

Darcy-Weisbach  friction  factor 

species  mass  fraction 

total  enthalpy,  J/kg 

specific  enthalpy,  J/kg 

heat  transfer  coefficient,  J/s/m2/K 

unit  vectors  for  the  cartesian  coordinates 

cell  index 

viscous  loss  coefficient 
coefficient  of  thermal  conductivity 
length,  m 
Mach  number 
molecular  weight 

mass  of  fluid  in  a  Lagrangian  cell,  kg 

piston  mass,  kg 

direction  cosine 

driver  gas  polytropic  index 

unit  vectors  for  the  cell  interface 

pressure,  Pa 

initial  partial  till  pressure  of  Helium  in  driver  tube 
initial  partial  till  pressure  of  Argon  in  driver  tube 
Prandtl  number 

primary  diaphragm  rupture  pressure,  MPa 

source  vector  in  the  gas-dynamic  equations 

(Lld2)  heat  transfer  rate,  J/s 

(mbcns2)  heat  flux,  W/m2 

gas  constant,  J /kg  •  K 

universal  gas  constant,  8.314  J/mole  •  K 

Reynolds  number 

tube  or  duct  radius,  m 

radial  coordinate,  m 

control  surface  of  the  cell 

Stanton  number 

Riemann  invariants 
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T  :  temperature,  K 

t  :  time,  s 

U  :  state  vector  in  the  gas-dynamic  equations 

u  :  velocity,  m/s 

Up  :  piston  velocity,  m/s 

urupt  '■  piston  velocity  at  diaphragm  rupture,  m/s 

Uref  :  reference  piston  speed  for  constant  driver  pressure 

V  :  cell  volume,  compression  tube  volume,  m3 

V  :  piston  velocity,  m/s 

v  :  volume,  m3 

w  :  work/unit  time  done  by  the  wall  shear  stress,  J/s 

X ,  x  :  piston  position,  m 

x,y,z  :  cartesian  coordinates,  m 

Z  :  intermediate  variable  for  the  Riemann  solver 

a  :  weighting  parameter 

f3  :  stretching  parameter 

f3  :  piston  over-driver  parameter  (=  urupt/Uref ) 

7  :  ratio  of  specific  heats 

A±  :  intermediate  variable  for  interpolation 

e  :  absolute  size  of  pipe  roughness  elements 

A  :  compression  ratio  for  the  free-piston  driver 

A  :  second  coefficients  of  viscosity,  Pa.s 

A  :  compressibility  factor 

fj,  :  viscosity,  Pa.s;  friction  coefficient 

7 r  :  3.14159... 

p  :  density,  kg/m3 

7  :  ratio  of  specific  heats 

t  :  wall  shear  stress,  Pa 

Q  :  recovery  factor 
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Subscripts 

A  :  reservoir 

A,  0  :  reservoir  initial  condition 

acc  :  acceleration  tube 

aw  :  adiabatic  wall  condition 

B  :  back  of  piston 

D  :  driver  tube 

D,  0  :  driver  tube  initial  condition 

exp  :  experimentally  measured  quantity 

F  :  front  of  piston 

/  :  friction  value 

i  :  cell  index,  inviscid 

j  :  cell  index 

j  ±  ^  :  interface  indices 

L,  R  :  left  and  right  states  for  the  Riemann  solver  or  flux  calculator 
loss  :  pipe  fitting  value 

m  :  quantity  at  piston  inflection  point  when  acceleration  is  zero 

max  :  maximum  value 

n  :  normal  to  the  cell  interface 

p  :  tangent  to  the  cell  interface 

:  piston 

s  :  nozzle  supply  (stagnation)  condition 

/  species  index 

sec  :  secondary  driver  tube 

shk  :  shock  tube 

tt  :  test  time 

v  :  viscous 

wall  :  wall  condition 

x,  y  :  coordinate  directions 

0  :  stagnation  property 

1,2  :  pre-  and  post-incident  shock  conditions  in  the  shock  tube 

Superscripts 

*  :  intermediate  state  for  the  Riemann  solver 

/  Eckert  reference  conditions 
(•  •  •)  :  cell  average 
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